A PRESSURE-BASED  COMPOSITE  GRID  METHOD  FOR 
COMPLEX  FLUID  FLOWS 


By 


JEFFREY  A.  WRIGHT 


A DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 
OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 
OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 
DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 


1993 


ACKNOWLEDGEMENTS 


I can  only  begin  to  express  my  sincere  gratitude  to  my  advisor,  Professor  Wei  Shyy, 
for  his  guidance,  never-ending  encouragement  and  unrelenting  humor  during  my  time 
working  with  him.  The  atmosphere  which  he  has  created  between  him  and  his  students  is  a 
special  one,  and  I truly  feel  that  he  is  not  only  a mentor,  but  a friend  as  well.  I also  wish  to 
thank  him  for  his  generous  financial  support,  which  allowed  me  to  maintain  a decent 
standard  of  living  during  the  preparation  of  this  work. 

I would  also  like  to  acknowledge  my  parents,  whose  unwavering  support, 
understanding  and  love  are  always  a source  of  encouragement,  and  whose  dedication  and 
foresight  gave  me  the  opportunity  to  pursue  this  path  in  the  first  place.  To  my  officemate 
Siddharth  (ST),  I can  only  say  that  I could  not  have  been  more  fortunate  to  have  such  a good 
friend,  and  I hope  time  and  space  cannot  separate  us.  I would  also  like  to  express  thanks  to 
my  good  friend  Uday,  who  is,  quite  simply,  one  of  the  finest  individuals  I have  ever  known. 
Thanks  are  also  due  to  ST’s  wife  Jyoti,  who  provided  many  a good  meal  and  taught  me  some 
of  the  finer  details  of  cooking  Indian  cuisine. 

I would  also  like  to  thank  the  members  of  my  committee.  Professors  David 
Mikolaitis,  Chen-Chi  Hsu,  Renwei  Mei,  Craig  Saltiel  and  especially  Professor  Richard 
Fearn,  who  opened  the  door  to  many  opportunities  for  me  in  the  beginning  of  my  research 
career.  Thanks  are  also  due  to  Madhukar  Rao  for  his  many,  and  often  useful, 
computer-related  suggestions  and  gadgets,  which  have  admittedly  saved  countless  hours  of 
my  time.  The  many  helpful  discussions  with  Jian  (James)  Liu,  concerning  composite  grid 
methods  are  also  appreciated. 


11 


TABLE  OF  CONTENTS 


Cage 

ACKNOWLEDGEMENTS  ii 

ABSTRACT  v 

CHAPTERS 

1 INTRODUCTION  1 

2 GOVERNING  EQUATIONS  AND  NUMERICAL  ALGORITHM  12 

2.1  Governing  Equations  and  Sequential  Solution  Procedure  12 

2.2  Discretization  of  the  Momentum  Equations 13 

2.3  Discretized  Form  of  the  Pressure  Correction  Equation  16 

3 COMPOSITE  GRID  METHOD  19 

3.1  Organization  of  Composite  Grids  19 

3.2  Global  Conservation  Conditions  for  the  Staggered  Grid  23 

3.3  Conservation  Strategy  for  Pressure-Based  Method  28 

3.4  Explicit  Conservation  Procedure  31 

3.5  Conditions  for  Local  and  Global  Conservation  37 

3.6  Numerical  Validation:  Lid-Driven  Cavity  Flow  41 

3.7  Computation  in  Multiply-Connected  Flow  Domains  46 

3.8  Implementaion  of  Higher-Order  Convection  Schemes  50 

3.9  Concluding  Remarks  55 

4 COMPOSITE  MULTIGRID  57 

4.1  Background  57 

4.2  Implementation  Issues  Associated  with  Composite  Multigrid  ^ 

4.3  Evaluation  of  Composite  Multigrid  Technique  72 

4.4  Concluding  Remarks  88 

iii 


5 COMPOSITE  GRID  METHOD  FOR  SCALAR  TRANSPORT  89 

5.1  Introduction  89 

5.2  Composite  Grid  Method  for  Natural  Convection  Flows 92 

5.3  Application  to  Double-Diffusive  Flow  in  a Solar  Pond 97 

5.4  Concluding  Remarks  126 

6 SUMMARY  AND  SUGGESTIONS  FOR  FUTURE  WORK  127 

REFERENCES  131 

BIOGRAPHICAL  SKETCH  135 


IV 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

A PRESSURE-BASED  COMPOSITE  GRID  METHOD  FOR 
COMPLEX  FLUID  FLOWS 

By 

Jeffrey  A.  Wright 
December  1993 


Chairman  : Wei  Shyy 

Major  Department : Aerospace  Engineering,  Mechanics  and  Engineering  Science 

In  this  work,  a pressure-based  composite  grid  algorithm  is  developed  for  solving  the 
incompressible  Navier-Stokes  equations  on  domains  composed  of  an  arbitrary  number  of 
overlapping  subgrids.  The  governing  equations  are  solved  independently  within  each  of  the 
subgrids  and  a staggered  grid  is  used  to  store  the  dependent  variables.  The  primary  emphasis 
is  on  the  development  of  a conservative  interface  treatment  for  transferring  information 
between  the  subgrids  of  the  system.  The  differences  between  the  conservative  interface 
scheme  developed  for  a staggered  grid  and  that  for  a nonstaggered  grid  are  discussed.  An 
organizational  scheme  is  developed  in  order  to  provide  a general  and  more  flexible  means 
for  handling  arbitrarily  overlain  subgrids.  In  order  to  provide  improved  computational 
efficiency,  a multigrid  technique  previously  implemented  for  a single  grid  solution 
algorithm  has  been  adopted  for  use  with  the  composite  grid  algorithm.  In  the  current 
procedure,  the  composite  grid  method  operates  as  an  independent  module  within  the 
framework  of  the  outer  multigrid  shell.  Results  for  simple  configurations  show  that 
composite  gridding  has  little  detrimental  effect  on  overall  multigrid  performance,  and  that 


even  for  geometrically  complex  flows  the  current  composite  multigrid  method  can  largely 
maintain  its  effectiveness.  With  the  extension  of  the  current  method  to  natural  convection 
flows  it  is  demonstrated  that  a locally  conservative  procedure  for  exchanging  information 
between  subgrids,  for  the  energy  and/or  species  equations,  can  allow  a non-unique  solution 
jump  across  block  interfaces.  A new  scheme,  based  on  a linear  interpolation  with  global 
correction  is  employed  which  allows  the  proper  solution  continuity  across  internal  block 
interfaces  to  be  maintained.  To  demonstrate  the  overall  effectiveness  of  the  technique,  the 
double-diffusive  flow  in  a solar  pond-type  configuration  is  computed.  With  an  initially 
gravitationally  stable  system  consisting  of  a two-component  linear  stratification  of  both 
temperature  and  salinity,  sidewall  heating  is  introduced,  resulting  in  an  instability  which 
causes  convective  cells  to  form  at  the  heated  wall  which  subsequently  propagate  into  the 
domain.  Two  classes  of  flows  resulting  from  different  initial  conditions  of  vertical  and  lateral 
system  stability  have  been  simulated.  Comparisons  are  made  with  existing  experimental 
simulations,  which  indicate  that  the  current  composite  grid  method  can  be  combined  with 
locally  refined  grids  to  provide  an  accurate  and  efficient  technique  for  tracking  complex, 
localized,  time-evolving  flows. 
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CHAPTER  1 
INTRODUCTION 


In  the  field  of  computational  fluid  dynamics  (CFD),  much  research,  within  the  past 
decade,  has  been  focused  toward  computing  complex  fluid  flows  of  practical  interest  to  the 
engineering  and  science  community.  Work  has  been  undertaken  in  several  key  areas  in  an 
attempt  to  understand  and  quantify  the  various  issues  which  must  be  resolved  before  real 
flows  can  be  handled  accurately  and  effectively  in  a generalized  manner.  Some  of  these  areas 
include  the  development  of  improved  models  for  describing  key  physical  mechanisms,  such 
as  turbulence,  combustion,  and  phase  change,  the  development  of  accurate  and  robust 
numerical  schemes  for  handling  convection  and  sharp  gradients  within  flowfields,  and  the 
implementation  of  numerical  algorithms  within  parallel  processing  environments.  One  area 
of  particular  concern,  however,  on  which  a large  brunt  of  the  load  hinges,  involves  the  ability 
to  handle  the  physical  geometric  complexity  commonly  found  in  many  engineering 
applications.  Since  the  ability  to  numerically  solve  systems  of  partial  differential  equations 
accurately  within  a geometrically  complex  domain  is  largely  dependent  on  the  quality  of  the 
grid  used  in  discretizing  the  domain,  the  success  of  the  overall  attempt  in  handling  real  flows 
largely  relies  on  the  ability  to  generate  suitable  grids  for  arbitrarily  complex  geometries. 

In  terms  of  handling  configurations  with  complex  geometry,  two  fundamentally 
different  approaches  have  been  developed  and  advocated  in  the  CFD  literature,  namely, 
unstructured  grids  and  composite  structured  grids.  With  unstructured  grids,  the  domain  is 
typically  discretized  with  a single  grid  composed  of  an  array  of  cells  which  can  have  an 
arbitrary  connectivity.  Usually  the  grid  is  composed  of  a single  cell  type,  typically  either 
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triangular  or  quadrilateral  (in  two  dimensions),  and  tetrahedral  (in  three  dimensions); 
however,  this  is  not  a requirement,  and  implementations  have  been  devised  which  use  a 
combination  of  different  cell  types  (Holmes  & Connell  1989).  Composite  structured  grid 
techniques,  on  the  other  hand,  as  the  name  indicates,  combine  a number  of  structured 
subgrids  to  form  the  complete  grid,  each  having  the  usual  i,j,k  ordering  of  cells.  For  many 
problems  of  engineering  interest,  the  generation  of  a single  structured  grid  which  adequately 
discretizes  the  domain  for  resolving  the  various  flow  features,  may  be  difficult  and  in  many 
cases  impractical.  The  composite  structured  grid  method  arose  out  of  the  need  for  handling 
flows  with  complex  geometry  while  maintaining  the  wealth  of  numerical  techniques  which 
have  been  developed  over  the  years  for  individual  structured  grids. 

Each  of  the  two  techniques  mentioned  above  has  its  own  advantages  and 
disadvantages.  Due  to  the  arbitrary  connectivity  of  the  cells,  unstructured  grids  provide  a 
degree  of  flexibility  in  discretizing  complex  geometries  which  cannot  be  obtained  with 
structured  grids.  With  the  development  of  automated  unstructured  grid  generation 
techniques  such  as  the  advancing  front  method  (Peraire  et  al.  1987,  Lohner  & Parikh  1988, 
Lohner  1988,  Morgan  et  al.  1991),  discretization  of  even  the  most  complex  domains  can  be 
handled  in  an  efficient  and  relatively  straightforward  manner.  In  addition  to  providing  a 
flexible  pre-solution  grid  generation  capability,  computational  algorithms  for  unstructured 
grids  are  also  ideally  suited  for  adaptive  grid  generation  techniques  involving  local 
enrichment.  Still,  however,  in  terms  of  grid  generation  alone,  several  issues,  such  as  the 
efficient  discretization  of  regions  containing  thin  layers,  such  as  boundary  layers,  and  more 
generally  the  ability  to  accurately  control  the  aspect  ratio  and  directionality  of  the  grid  cells, 
need  to  be  resolved  before  unstructured  grids  can  be  viewed  as  completely  adequate  for 
arbitrary  configurations  requiring  the  solution  of  the  Navier-Stokes  equations.  In  this  regard, 
techniques  employing  either  a combined  unstructured/structured  grid  approach  or  a mixed 
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unstructured  element  approach  (Holmes  & Connell  1989)  have  been  developed  which  can 
more  efficiently  accommodate  flow  regions  (most  notably  boundary  layers),  where  a 
directional  length  scale  disparity  is  imposed. 

As  with  many  numerical  techniques  where  some  beneficial  property  is  gained,  some 
sacrifice  must  be  made.  Due  to  the  lack  of  inherent  order  of  the  unstructured  grid, 
significantly  more  operations  per  grid  node  are  required  compared  with  structured  grids  to 
assemble  the  linear  system  of  equations  which  form  the  discrete  analogue  of  the  governing 
system  of  equations  (Lohner  1987).  This  is  primarily  due  to  the  required  gathering/scattering 
operations  which  involve  many  additional  indexing  operations,  which  are  not  required  for 
structured  grids.  Furthermore,  additional  data  structures  are  required  to  store  the 
cell-to-node  and  face-to-node  connectivity  information  which  is  implicitly  contained  within 
structured  grids.  Therefore,  with  unstructured  grids,  a sacrifice  is  felt  in  both  the  overall 
numerical  efficiency  achieved  in  obtaining  solutions  and  in  the  use  of  computer  memory 
resources.  For  a complete  review  of  the  relative  merits  and  demerits  of  unstructured  grid 
techniques  and  the  state-of-the-art  methods  now  employed  in  the  CFD  community,  the 
reader  is  referred  to  Lohner  (1987)  and  Zienkiewicz  (1992). 

In  general,  composite  structured  grids  are  less  flexible  than  unstructured  grids  in 
resolving  highly  complex  geometries;  however,  for  many  configurations  of  interest, 
composite  structured  grids  prove  sufficient  for  providing  an  accurate  discretization  of  the 
domain  (see,  for  example,  Stewart,  1991).  For  composite  structured  grid  algorithms  the 
difficulties  incurred  in  the  grid  generation  process,  for  most  configurations,  are  more  than 
outweighed  by  the  overall  advantage  which  is  obtained  in  terms  of  both  numerical  accuracy 
and  efficiency  as  well  as  the  efficient  utilization  of  computer  memory  resources.  Additional 
complexity  is,  however,  introduced  due  to  the  presence  of  the  artificial  boundaries  separating 
the  various  subgrids  (blocks)  which  are  used  to  construct  the  complete  domain.  The  level 
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of  additional  complexity  is  also  highly  dependent  on  the  specific  composite  grid  technique 
used,  as  well  as  the  method  employed  for  updating  the  boundary  information  for  the  various 
subgrids.  Over  the  years,  several  distinct  composite  grid  techniques  have  evolved, 
employing  either  patched  (abutting)  grids  or  overlaid  grids  (which  include  the  so-called 
Chimera- type  grids).  In  general,  techniques  employing  overlaid  grids  provide  an  advantage 
over  those  which  use  patched  grids,  in  terms  of  the  ability  to  discretize  arbitrarily  complex 
domains;  however,  the  price  paid  for  this  additional  flexibility  is  usually  felt  in  the 
implementation  difficulties  which  are  incurred  in  obtaining  internal  boundary  information 
for  the  more  generally  oriented  subgrids.  For  a concise  review  of  the  various  composite  grid 
techniques  which  are  in  use,  see  Steger  and  Benek  (1987). 

For  any  composite  grid  algorithm  in  which  the  governing  equations  are  solved 
independently  within  each  of  the  subgrids,  one  of  the  primary  issues  of  importance  concerns 
the  transfer  of  information  between  the  various  subgrids  during  the  course  of  the  outer 
iteration  process  in  which  one  sweeps  sequentially  from  one  subgrid  to  the  next.  Strategies 
for  transferring  information  between  subgrids  can  be  generally  classified  as  either 
conservative  or  nonconservative.  Nonconservative  transfer  schemes  usually  employ  some 
form  of  direct  interpolation  in  order  to  obtain  the  internal  boundary  values  for  each  of  the 
subgrids  (Atta  & Vadyak  1983,  Chesshire  & Henshaw  1990,  Steger  1991).  These  schemes 
are  the  most  primitive  type  of  transfer  scheme  and  are  usually  employed  for  composite  grid 
techniques  in  which  arbitrarily  complex  subgrid  overlappings  are  allowed  (such  as  the 
Chimera  technique),  for  no  other  reason  than  the  fact  that  the  implementation  of  conservative 
schemes  is  more  difficult.  Such  transfer  schemes  are  sufficient,  in  terms  of  just  obtaining 
solutions,  for  finite  difference  techniques  which  do  not  correspond  to  a particular  finite 
volume  formulation  of  the  governing  equations;  however,  they  can  preclude  converged 
solutions  from  being  obtained  with  finite  volume  solvers  unless  an  extra  degree  of  freedom 
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is  provided  by  the  boundary  conditions  for  ensuring  global  conservation,  or  global 
conservation  is  explicitly  enforced  in  some  manner. 

For  the  application  of  finite  volume  techniques  to  general  composite  grid  systems, 
some  form  of  conservative  internal  boundary  treatment  should  be  adopted.  Conservative 
transfer  schemes  can  be  further  classified  according  to  both  their  local  and  global 
conservation  properties.  A number  of  schemes  have  been  developed  which  exhibit  strict 
global  conservation,  but  not  local  conservation.  Most  of  these  schemes  are  based  on  an 
interpolation  procedure  combined  with  a global  correction  strategy  (Davis  & Thompson 
1992,  Tu  & Fuchs  1992).  Due  to  the  use  of  interpolation,  any  concept  of  local  conservation 
is  destroyed;  however,  global  conservation  can  still  be  maintained  by  computing  the  net 
global  imbalance  of  the  appropriate  conserved  quantity  across  the  internal  boundaries  of  a 
typical  subgrid  and  then  redistributing  an  equal  and  opposite  amount  across  these  boundaries 
in  some  manner.  In  order  to  maintain  accuracy  with  this  strategy  in  situations  in  which 
isolated  features  containing  sharp  gradients  pass  through  the  internal  boundaries,  effort  must 
be  made  to  ensure  that  the  redistribution  is  performed  in  an  equitable  manner,  whereby 
regions  of  the  interface  which  incur  larger  conservation  errors  due  to  interpolation  receive 
a larger  proportion  of  the  flux  which  is  redistributed  to  the  internal  boundary  in  the  global 
correction  step  (Davis  & Thompson  1992). 

Fully  conservative  interface  schemes,  exhibiting  both  local  and  global  conservation 
properties,  have  also  been  developed,  mainly  in  the  context  of  density-based  methods  for 
compressible  fluid  flows  (Rai  1986a,  1986b).  With  these  schemes,  local  conservation  is 
strictly  enforced  and  is  used  as  the  basis  for  also  ensuring  global  conservation.  Local 
conservation  along  an  internal  boundary  is  ensured  by  directly  computing  boundary  fluxes 
for  each  of  the  control  volumes  lying  along  the  boundary  in  a manner  which  is  consistent 
with  the  representation  of  these  fluxes  within  the  neighboring  subgrids  which  provide  the 
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flux  information.  The  use  of  local  conservation  as  the  framework  for  providing  internal 
boundary  conditions  ensures  that  no  redistribution  scheme  is  required,  and  thus  the  errors 
associated  with  interpolation/correction  schemes  do  not  appear  with  fully  conservative 
schemes.  Results  have  been  presented  for  compressible  flow  simulations  which  demonstrate 
that  quite  accurate  results  can  be  obtained  with  fully  conservative  techniques,  even  for  cases 
where  complex  fluid  flow  features  pass  directly  through  the  internal  boundaries  of  the  grid 
system  (Rai  et  al.  1984). 

While  substantial  work  has  been  conducted  in  developing  fully  conservative 
interface  schemes  for  density-based  algorithms  applied  to  compressible  flows,  relatively 
little  attention  has  been  paid  in  developing  analogous  schemes  for  incompressible  flows 
where  pressure-based  algorithms  are  commonly  employed.  Several  works  have  appeared  in 
the  literature;  however,  so  far,  the  composite  grid  schemes  which  have  been  developed  have 
employed  either  direct  interpolation  (Meakin  & Street  1988)  or  interpolation  with 
conservative  correction  (Tu  & Fuchs  1992).  In  a recent  work  by  Perng  and  Street  (1991), 
a composite  grid  technique  was  developed  for  the  incompressible  Navier-Stokes  equations 
employing  a pressure-based  method  with  a staggered  grid  arrangement;  however,  the 
interface  treatment  developed  there  employs  artificial  pressure  boundary  conditions  (of  both 
the  Dirichlet  and  Neumann  type)  at  the  internal  boundaries  and  is  not  conservative.  Thus, 
it  appears  that  no  concerted  effort  has  been  made  in  exploring  the  issues  involved  with  the 
development  of  fully  conservative  interface  schemes  for  pressure-based  algorithms  applied 
to  incompressible  flows. 

In  this  work,  a composite  grid  procedure  is  developed  for  solving  the  steady, 
incompressible  Navier-Stokes  equations  on  domains  comprised  of  an  arbitrary  number  of 
overlapping  subgrids.  A pressure-correction  algorithm  in  the  spirit  of  SIMPLE  (Patankar 
1980)  is  used  in  conjunction  with  a staggered  grid  system  to  solve  the  continuity  and 
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momentum  equations  in  a sequential  fashion  within  each  of  the  subgrids  composing  the 
domain  and  a fully  conservative  interface  treatment  is  devised  for  transferring  information 
between  the  various  subgrids.  There  are  a number  of  merits  gained  from  using  a sequential 
solution  technique  like  SIMPLE  along  with  a staggered  grid.  Sequential  solution  techniques 
allow  one  to  accommodate  a different  number  of  equations  depending  on  the  physics 
involved,  without  the  need  for  reformulating  the  solution  algorithm.  In  addition,  the  basic 
algorithm  can  be  extended  in  a unified  framework  to  handle  the  various  flow  regimes,  from 
incompressible  to  hypersonic  (Shyy  & Braaten  1988,  Rhie  1989).  Some  of  the  merits  of  the 
staggered  grid  include  the  compactness  of  the  discretization  formulae,  the  elimination  of  the 
need  for  higher-order  smoothing  to  prevent  odd-even  decoupling  of  the  pressure  field,  as 
well  the  elimination  of  the  need  for  artificial  pressure  boundary  conditions.  While  the  use 
of  the  staggered  grid  increases  the  complexity  of  the  numerical  scheme,  in  the  author’s 
opinion,  the  basic  merits  of  the  staggered  grid,  mentioned  above,  outweigh  the  drawbacks. 

In  the  following,  a brief  review  of  the  governing  equations  is  first  presented,  along 
with  a summary  of  the  sequential  solution  process  which  is  used  to  solve  the  discretized  form 
of  these  equations  on  single  structured  grid  domains.  Next,  an  in-depth  discussion  is  given 
regarding  the  development  of  the  current  composite  grid  approach  for  the  pressure-based 
algorithm.  The  first  issue  of  importance  in  this  area  concerns  the  organization  of  composite 
grids  comprised  of  an  arbitrary  number  of  overlapping  subgrids.  Such  systems  can  become 
very  complex,  especially  in  regions  where  multiple  subgrids  overlap,  and  therefore,  some 
strategy  must  be  developed  to  consistently  organize  these  systems  so  that  a conservative 
internal  boundary  scheme  can  be  implemented  in  a straightforward  manner.  Once  an 
organizational  scheme  has  been  proposed,  the  global  conservation  conditions  for  the 
staggered  grid  are  detailed,  followed  by  a discussion  of  the  conservative  strategy  for  the 
pressure-based  method.  Finally,  a detailed  description  of  the  explicit  conservation  procedure 
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used  for  maintaining  local  conservation  is  presented  and  the  overlapping  constraints 
required  for  simultaneously  ensuring  both  local  and  global  conservation  are  discussed. 

Following  the  development  of  the  composite  grid  technique,  several  model  problems 
are  investigated  in  order  to  both  validate  the  computational  procedure,  and  demonstrate  the 
capability  of  the  organizational  scheme  in  handling  configurations  exhibiting  marked 
geometric  complexity.  The  test  problem  used  for  numerical  validation  is  the  well-known 
square,  lid-driven  cavity  flow.  For  this  problem,  a composite  grid  case  is  computed  to 
numerically  verify  the  computational  algorithm,  as  well  as  to  demonstrate  the  computational 
efficiency  which  can  be  obtained  with  the  composite  grid  technique,  in  handling  internal 
flows  exhibiting  complex  recirculation  regions.  The  next  case  is  representative  of  a typical 
geometrically  complex,  multiply-connected  flow  domain  and  consists  of  a channel  with 
internal  baffles.  This  configuration  was  originally  inspired  by  an  electronics  application, 
where  the  internal  baffles  represent  electronic  components  housed  in  an  outer  cabinet 
through  which  a cooling  fluid  is  being  pumped.  The  computed  results  for  this  flow  confirm 
that  the  current  composite  grid  procedure  is  indeed  capable  of  resolving  complex  geometries 
similar  to  those  found  in  many  real-world  applications  of  practical  interest. 

For  any  computational  algorithm  employing  relaxation  techniques  to  solve  the  linear 
system  of  equations  resulting  from  the  discretization  of  the  governing  flow  equations,  the 
use  of  multigrid  techniques,  whereby  the  discrete  equations  are  solved  on  a sequence  of  grids 
of  varying  resolution,  can  substantially  increase  the  rate  of  convergence  of  the  iterative 
solution  process  (by  an  order  of  magnitude  or  more  in  some  cases).  The  multigrid  idea, 
originally  proposed  by  Brandt  (1977),  has  been  successfully  applied  within  the  area  of  CFD 
in  recent  years  for  a wide  variety  of  fluid  flow  computations,  including  the  solution  of  highly 
nonlinear  systems  of  equations  such  as  the  Navier-Stokes  equations.  Both  single  grid  (Ghia 
et  al.  1982,  Thompson  & Ferziger  1989,  Bruneau  & Jouron  1990)  and  composite  grid 
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(Hinatsu  & Ferziger  1991,  Tu  & Fuchs  1992)  solution  frameworks  have  been  devised  and 
tested,  with  varying  degrees  of  success  reported,  depending  on  the  governing  equations 
being  solved  as  well  as  the  particular  multigrid  implementation  employed.  In  this  work,  a 
multigrid  technique  previously  employed  by  Shyy  and  Sun  (1993)  for  solving  the 
Navier-Stokes  equations  within  single-grid  domains  has  been  adopted  and  implemented 
within  the  framework  of  the  current  composite  grid  procedure. 

In  extending  the  multigrid  technique  to  composite  grids,  several  key  implementation 
issues  arise,  which  are  critical  to  the  success  of  the  overall  procedure.  One  of  these  issues 
concerns  the  procedural  relationship  between  the  composite  grid  procedure  and  the 
multigrid  strategy.  Several  distinct  implementations  are  possible,  each  having  a unique 
coupling  relationship  between  composite  grid  and  multigrid.  In  order  to  provide  a robust 
capability  for  handling  arbitrarily  complex  flows,  a methodology  is  established  here, 
whereby  the  composite  grid  procedure  serves  as  an  independent  module  operating  within  the 
outer  framework  of  the  multigrid  cycling  algorithm.  Another  issue  of  importance  revolves 
around  the  choice  of  the  grid  coarsening  strategy  used  within  the  multigrid  cycling  step.  For 
composite  grid  techniques  employing  overlaid  grids,  certain  minimum  overlapping 
requirements  must  be  observed  in  order  for  internal  boundary  information  to  be  successfully 
communicated  among  the  various  subgrids  of  the  system.  These  requirements  lead  to  direct 
implications  regarding  the  type  of  coarsening  strategy  which  is  most  suitable  for  a combined 
composite  multigrid  procedure.  These  implications  are  discussed  in  detail,  and  from  the 
requirements  imposed  by  the  composite  procedure,  a general  multigrid  coarsening  strategy 
is  adopted  which  ensures  that  minimum  overlapping  requirements  are  met  at  all  multigrid 
levels. 

Following  the  development  of  the  basic  composite  multigrid  method,  a number  of 
test  cases  are  presented  to  determine  the  effect  of  composite  gridding  on  overall  multigrid 
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performance.  Composite  grid  test  cases  for  the  square,  lid-driven  cavity  are  first  used  here 
to  compare  the  relative  performance  of  composite  multigrid  to  that  obtained  with 
conventional  single-grid  implementations  for  simple  configurations.  To  assess  the 
robustness  of  the  current  capability  in  handling  geometrically  complex  flows,  the  composite 
multigrid  technique  is  applied  to  two  multiply-connected  flow  domains.  Results  for  these 
test  cases  indicate  that  the  current  technique  can  largely  maintain  its  efficiency  even  for 
geometrically  complex  recirculating  flows. 

In  the  last  section  of  this  work,  some  issues  related  to  the  solution  of  natural 
convection  flows  with  composite  grid  methods  are  examined.  The  governing  equations  of 
concern  are  the  time-dependent  incompressible  Navier-Stokes  equations  with  the 
Boussinesq  approximation  employed,  along  with  the  corresponding  transport  equations  for 
additional  scalar  quantities  affecting  the  density  (such  as  temperature  and  salinity).  The 
momentum  equations  are  identical  to  those  employed  in  the  previous  sections  with  the 
exception  of  the  gravitational  force  term  in  the  v-momentum  equation.  Since  this  is  a body 
force  term,  it  has  no  effect  on  the  conservation  procedure  previously  developed,  and  thus, 
the  only  additional  requirement  at  the  internal  boundaries  is  to  develop  an  information 
transfer  strategy  for  the  auxiliary  equations  for  temperature,  salinity,  etc.  In  this  regard,  it 
is  demonstrated  that  a locally  conservative  flux  estimation  based  strictly  upon  values 
obtained  within  neighboring  grid  blocks  allows  an  arbitrary  jump  across  block  boundaries 
which  is  a function  of  the  local  flow  conditions.  The  possibility  of  the  existence  of  such  a 
jump  is  demonstrated  using  a one-dimensional  model  problem  and  is  also  confirmed 
numerically  for  a two-dimensional  natural  convection  cavity  flow.  Based  upon  this  finding, 
a new  boundary  scheme  based  on  a linear  interpolation  with  correction  procedure  is  devised 
which  allows  the  proper  solution  continuity  to  be  maintained  across  internal  boundaries  for 


natural  convection  flows. 
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As  a application  of  the  overall  composite  grid  method  for  natural  convection  flows, 
a time-dependent  simulation  of  a double-diffusive  flow  in  a rectangular  enclosure  is 
investigated.  With  an  initially  gravitationally  stable  system  consisting  of  a two-component 
linear  stratification  of  both  temperature  and  salinity,  sidewall  heating  is  introduced,  resulting 
in  an  instability,  which  causes  concentrated  convective  cells  to  form  at  the  heated  wall  and 
then  propagate  into  the  domain.  Study  of  this  problem  was  initially  inspired  by  observations 
of  well-defined  stable  convective  layers  in  natural  reservoirs  such  as  lakes  and  oceans  (Hoare 
1966,  Huppert  & Turner  1980).  Similar  observations  have  also  been  made  during  the 
operation  of  solar  energy  collection  salt  ponds,  where  sidewall  heating  due  to  the  angle  of 
incident  radiation  of  the  sun  was  thought  to  induce  horizontal  convection  patterns  which 
eventually  destroyed  the  stabilizing  salt  stratification  of  the  pond,  and  thus  its  solar  energy 
collection  capability  (Akbarzadeh  & Manins  1988,  Sherman  & Imberger  1991). 

Another  point  of  interest  to  be  demonstrated  here  is  the  effectiveness  with  which  the 
current  composite  grid  method  can  be  utilized  in  computing  time-dependent  flows 
containing  localized  flow  features  with  sharp  gradients.  Here,  an  adaptive  composite  grid 
technique  based  on  local  grid  refinement  is  employed  to  track  the  sidewall  intrusions  as  they 
develop  and  propagate  into  the  domain.  The  initial  grid  consists  of  a fine  grid  in  the  initiation 
region  near  the  heated  wall  and  a coarser  grid  throughout  the  rest  of  the  domain.  As  the 
intrusions  develop  and  propagate  away  from  the  wall,  the  fine  grid  expands  into  the  domain 
to  match  the  speed  of  the  intrusion  front,  while  the  coarser  grid  is  simultaneously  retracted 
at  the  same  speed.  In  this  case,  the  coarse  grid  is  acting  essentially  as  a connection  to  the 
far-field  boundary  condition.  Comparisons  with  experimental  results  indicate  that  the 
current  method  can  provide  an  accurate  and  efficient  technique  for  tracking  localized, 
time-evolving  features.  Finally,  a summary  of  the  major  conclusions  of  the  work  is 
presented,  and  some  suggestions  for  future  work  in  this  area  are  given. 


CHAPTER  2 

GOVERNING  EQUATIONS  AND  NUMERICAL  ALGORITHM 


2.1  Governing  Equations  and  Sequential  Solution  Algorithm 


In  order  to  provide  a basis  for  the  development  and  discussion  of  the  composite  grid 
procedure  to  follow,  a brief  review  of  the  governing  equations,  sequential  solution 
technique,  and  discretization  procedure  is  presented  here.  For  additional  details,  the  reader 
is  referred  to  the  work  of  Patankar  (1980),  where  a complete  development  and  discussion 
of  the  SIMPLE  algorithm  is  given  for  single  grid  domains.  The  governing  equations  used 
here  are  the  steady,  incompressible,  two-dimensional  Navier-Stokes  equations  along  with 
the  corresponding  form  of  the  continuity  equation,  written  in  conservative  form  as  follows: 


According  to  SIMPLE,  an  initial  guess  is  first  made  for  the  velocity  and  pressure 
field.  The  momentum  equations  are  then  solved  to  yield  updated  estimates  for  the  velocity 
field.  Next,  a pressure  correction  equation  is  solved  to  give  an  updated  estimate  of  the 
pressure  field.  The  velocity  field  obtained  above  is  also  corrected  using  the  pressure 
corrections  to  obtain  a continuity  satisfying  velocity  field  at  each  iteration.  If  the  continuity 
and  momentum  equations  are  not  satisfied  to  the  desired  tolerance,  then  the  process  is 


(2.1a) 


(2.1b) 


(2.1c) 
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repeated  using  the  current  estimates  of  the  velocity  and  pressure  fields.  In  this  work,  a 
staggered  grid  system  is  adopted,  as  detailed  in  the  introduction.  The  particular  staggered 
grid  arrangement  chosen  here  is  shown  in  Figure  2.1.  In  the  figure,  the  physical  locations 
of  the  two  Cartesian  velocity  components,  pressure,  and  density  associated  with  the  grid 
point  ( i , j ) are  given.  With  the  staggered  grid  convention  adopted  in  Figure  2.1,  the 
unknowns  for  a typical  grid  appear  as  shown  in  Figure  2.2(a).  The  physical  locations  of  the 
known  boundary  velocity  components  are  shown  in  Figure  2.2(b).  No  pressure  boundary 
conditions  are  required  for  the  staggered  grid  arrangement. 

2.2  Discretization  of  the  Momentum  Equations 

Discretization  of  the  momentum  equations  is  accomplished  using  a control  volume 
approach.  Figure  2.3  shows  a finite  volume  representation  for  a typical  control  volume  about 
the  unknown  u velocity  component  up . Neighboring  unknowns  have  also  been  shown  in  the 
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Figure  2.2  Location  of  unknown  and  known  variables  for  the  staggered  grid, 
(a)  Unknown  velocity  components  and  pressures;  (b)  Known 
boundary  velocity  components. 

figure.  The  notation  shown  in  the  figure,  and  adopted  here  throughout,  is  standard  for  the 
SIMPLE  algorithm  with  the  staggered  grid  system.  Upon  integrating  the  u-momentum 
equation  over  the  u control  volume  limits,  one  arrives  at  the  general  form  of  the  conservation 
equation  for  u-momentum  as  follows: 

Apup  = ^ A,-  ut  + ipw  ~ PE ) Ay  (2.2) 

i = E,WJJ,S 

The  coefficients  AP  and  A,-  represent  the  combined  convection  and  diffusion  fluxes 
associated  with  each  of  the  nodes.  In  order  to  aid  the  discussion  later  needed  for  the 
development  of  the  composite  grid  algorithm,  representative  expressions  of  these  A;  terms 
are  given  as  follows,  where  the  second-order  central  difference  scheme  is  adopted  for  both 


the  convection  and  diffusion  terms. 
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Figure  2.3.  U control  volume  used  in  deriving  the  discretized  form  of  the 
u-momentum  equation. 
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Equations  (2.3)  are  part  of  a generalized  form  capable  of  representing  several  different 
discretization  schemes  (Patankar  1980,  Shyy  et  al.  1992,  Thakur  & Shyy  1993).  The  F’s  and 
D’s  represent  the  convection  and  diffusion  flux  coefficients  at  the  faces  of  the  control  volume 
and  are  given  by  the  following  expressions: 


In  the  above  expressions,  [[a  , bj  represents  the  maximum  of  the  two  arguments  a and  b. 
It  is  noted  that  the  terms  appearing  within  the  parentheses  of  Equation  (2.3e)  are  collectively 
equal  to  zero  because  of  the  mass  continuity  constraint.  They  are  explicitly  retained  for 
helping  the  understanding  of  the  interface  treatment  to  be  presented  later.  A similar 
procedure  can  be  performed  for  the  v-momentum  equation  by  integrating  the  v-momentum 
equation  over  a v control  volume,  resulting  in  the  following: 


where  the  coefficients  are  identical  to  those  previously  given  for  the  u-momentum  equation. 


The  continuity  and  momentum  equations  can  be  used  to  formulate  a pressure 
correction  equation.  The  pressure  correction  is  used  to  update  the  pressure  field  and  in 
conjunction  with  the  velocity  correction  formulas  to  obtain  a continuity  satisfying  velocity 
field  at  each  iteration  step.  In  order  to  aid  the  discussion  of  the  composite  grid  procedure  to 
be  developed  later,  the  discretized  form  of  the  pressure  correction  equation  is  presented 


Fe  = (, gu)edy  Fw  = (gu)  wAy  Fn  = (gv)nAx  Fs  = (pv)ydx  (2.4) 


(2.5) 


(2.6) 


i = E.W.N.S 


2.3  Discretized  Form  of  the  Pressure  Correction  Equation 


below. 


17 


app'p  = aEp'E  + awp'w  + aNp'N  + asp's  + b 

QeU ly)2  Qw(4y)2  QniAx)2  _ QS{Z. \x)2 

°E  ~ ~A  “w  ~ ~A  ~ ~A  as  ~ ~A 

APe  APw  APn  APs 

Clp  = dE  + dyy  + dpj  + dp 

b = [ ( QU*)W  - 0 ou)e  ]Ay  + [ (gv*)s  - (qv*)„  ]Ax 

In  the  above  expressions,  the  starred  velocity  components  represent  those  values  obtained 
from  the  most  recent  solution  of  the  momentum  equations.  The  A P terms  are  the  coefficients 
AP  in  the  discretized  forms  of  the  momentum  equations  for  the  velocity  components  located 
on  the  faces  of  the  pressure  prime  control  volume. 

It  should  be  noted  here  that  the  pressure  correction  equation  presented  above  is  not 
an  exact  representation  of  the  finite  volume  form  of  the  combined  continuity  and  momentum 
equations.  Truncated  forms  of  the  equations  relating  the  velocity  and  pressure  corrections 
are  used  in  the  derivation  of  Equation  (2.7)  in  order  to  obtain  a five-point  stencil.  In  Equation 
(2.7),  the  most  important  term  is  the  source  term  b.  This  source  term  represents  the  net 
imbalance  of  mass  flux  into  a pressure-correction  control  volume  and  is  used  to  drive  the 
equation  into  producing  pressure  corrections  which,  when  combined  with  the  velocity 
correction  formulas,  produce  a continuity  satisfying  velocity  field  at  the  end  of  each  outer 
iteration,  an  outer  iteration  being  defined  as  a sweep  through  the  u-momentum, 
v-momentum,  and  pressure-correction  equations.  In  this  regard,  many  different 
pressure-correction  equations  of  the  type  shown  in  Equation  (2.7)  can  be  formulated,  in 
terms  of  the  coefficients  a(- ; however,  any  of  these  formulations  must  satisfy  the  relation 
shown  in  Equation  (2.8b)  and  also  carry  a source  term  identical  to  that  shown  in  Equation 


(2.7) 

(2.8a) 

(2.8b) 

(2.9) 
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(2.9),  since  these  are  the  only  representations  which  will  produce  a uniform 
pressure-correction  field  when  mass  continuity  is  satisfied  for  all  control  volumes. 


CHAPTER  3 

COMPOSITE  GRID  METHOD 
3.1  Organization  of  Composite  Grids 

One  of  the  primary  difficulties  in  dealing  with  composite  grids  is  the  organizational 
task  of  determining  the  information  flow  from  block  to  block.  For  a composite  grid 
consisting  of  only  two  overlapping  grid  blocks,  or  for  any  composite  grid  in  which  no  more 
than  two  blocks  in  the  system  overlap  at  any  single  point  in  the  domain,  the  channel  of 
information  flow  is  already  determined,  since  there  is  only  one  neighboring  block  for  each 
block  in  the  system  which  can  provide  the  required  internal  boundary  data.  However,  for 
cases  in  which  three  or  more  grid  blocks  in  the  composite  grid  system  overlap,  determination 
of  the  internal  boundary  information  flow  paths  becomes  more  difficult.  In  these  cases, 
several  grid  blocks  may  be  available  to  provide  information  across  some  parts  of  an  internal 
boundary,  while  for  other  parts  of  the  boundary,  information  may  only  be  obtained  from  one 
other  neighboring  block.  An  example  of  this  may  be  seen  from  the  simple  composite  grid 
system  shown  in  Figure  3.1.  Here,  for  the  lower  side  of  block  two,  which  is  entirely  an 
internal  boundary,  for  the  right  portion  of  this  boundary,  information  can  be  obtained  only 
from  block  three;  however,  for  the  left  portion  of  this  boundary,  information  may  be  obtained 
from  either  block  one  or  block  three. 

The  problem  of  handling  composite  grids  consisting  of  multiple  overlapping  blocks, 
at  least  in  terms  of  information  transfer,  is  equivalent  to  determining  what  information  flows 
across  particular  segments  of  each  side  of  each  block  in  the  system.  Each  block  side  can  have 
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two  types  of  segments  associated  with  it,  which  will  be  designated  as  boundary  condition 
segments  and  internal  boundary  segments.  Boundary  condition  segments,  which  include 
both  Dirichlet  and  outflow  (i.e.  gradient  condition)  segments,  are  specified  at  the  outset  of 
the  problem,  and  in  terms  of  information  flow,  all  required  information  is  specified  via  the 
particular  boundary  condition  for  each  segment.  Internal  boundary  segments,  on  the  other 
hand,  are  created  by  the  overlappings  of  the  blocks  in  the  composite  grid.  Each  side  of  each 
block  in  the  system  is  spanned  by  a combination  of  internal  boundary  segments  and/or 
boundary  condition  segments.  Across  each  of  the  internal  boundary  segments,  information 
is  obtained  from  a single  neighboring  grid  block.  While  the  internal  boundary  segments  are 
essentially  determined  by  the  arrangement  of  the  composite  grid  system,  for  composite  grids 
containing  regions  in  which  more  than  two  blocks  overlap,  the  set  of  internal  boundary 
segments  is  not  unique.  This  again  relates  to  the  fact  that  information  may  be  obtained  across 
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portions  of  the  internal  boundaries,  from  more  than  one  of  the  neighboring  blocks.  Thus, 
some  strategy  must  also  be  developed  for  selecting  a unique  set  of  internal  boundary 
segments  for  each  of  the  blocks.  Once  the  internal  boundary  segments  for  each  block  in  the 
composite  grid  system  have  been  determined,  then  along  with  the  boundary  condition 
segments  of  the  blocks,  the  complete  paths  of  information  flow  into  each  block  are  specified, 
and  the  conservative  internal  boundary  scheme  can  be  implemented. 

The  procedure  used  for  determining  a unique  set  of  internal  boundary  segments  for 
each  block  in  the  composite  grid  system  can  be  summarized  as  follows.  First  an  intersection 
test  is  performed  for  each  block  in  the  composite  grid  with  every  other  block  to  arrive  at  a 
preliminary  set  of  internal  boundary  segments.  For  general  composite  grids,  the  internal 
boundary  segments  may  overlap  along  portions  of  the  boundaries,  resulting  in  an  ambiguity 
in  terms  of  information  transfer  there.  In  order  to  resolve  this  ambiguity,  each  block  in  the 
composite  grid  system  is  given  an  overlapping  priority,  so  that  when  two  internal  boundary 
segments  overlap,  the  segment  with  the  highest  priority  takes  precedence.  In  this  way,  a 
unique  set  of  internal  boundary  segments  can  be  obtained.  This  two-step  process  is  best 
illustrated  by  example. 

Consider  again,  the  composite  grid  system  shown  in  Figure  3.1.  We  will  determine 
a unique  set  of  internal  boundary  segments  for  block  two  of  this  three  block  composite  grid. 
First,  an  intersection  of  block  two  independently  with  each  of  blocks  one  and  three  is 
performed.  The  preliminary  internal  boundary  segments  obtained  from  these  intersection 
tests  are  shown  in  Figure  3.2(a).  The  number  by  each  of  the  segments  indicates  that  that 
segment  was  created  by  the  intersection  of  block  two  with  the  block  of  that  number.  For  the 
lower  internal  boundary  as  well  as  the  left  internal  boundary  of  block  two,  the  internal 
boundary  segments  are  seen  to  overlap  at  certain  locations.  Now,  the  following  overlapping 
priorities  are  further  specified  for  the  blocks  in  the  composite  grid.  Block  number  one  is 
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(a) 


Figure  3.2.  Determination  of  internal  boundary  segments  for  block  two 
of  composite  grid  of  Figure  3.1.  (a)  Preliminary  internal 
boundary  segments;  (b)  Final  internal  boundary  segments. 


given  the  highest  priority,  followed  by  block  two,  and  then  block  three.  After  using  this 
priority  information  to  eliminate  internal  boundary  segment  overlaps,  the  final  set  of  internal 
boundary  segments  appears  as  shown  in  Figure  3.2(b).  With  this  set  of  internal  boundary 
segments,  the  information  channels  through  the  internal  boundaries  are  precisely  specified. 
Along  with  the  boundary  condition  segments,  a complete  specification  of  the  information 
flow  into  block  two  is  obtained.  It  should  be  noted  that  this  procedure  for  determining  the 
exchange  of  information  within  a composite  grid  arrangement  is  quite  general  and  can  be 
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applied  to  composite  grids  created  from  any  number  of  grid  blocks  overlaid  in  an  arbitrary 
fashion. 


3.2  Global  Conservation  Conditions  for  the  Staggered  Grid 

Consider  the  u-momentum  equation  written  in  conservative  form  as 
follows: 

Ex  + Fy  = 0 (3.1a) 

where  the  terms  E and  F represent  the  local  fluxes  of  u momentum  in  the  x and  y 
directions  respectively  and  are  written  as 

E = Quu  - n + p F = guv  - (3.1b) 

Integrating  this  equation  over  a typical  u control  volume  yields 

(Ee  - Ew)Ay  + (Fn  - Fs)Ax  = 0 (3.2) 

The  terms  E^Ay  and  E^Ay  are  respectively  the  total  fluxes  of  u-momentum  through  the 
east  and  west  control  volume  faces  and  the  terms  F„Ax  and  F Ax  are  respectively  the  total 
fluxes  of  u-momentum  through  the  north  and  south  control  volume  faces. 

Now  consider  the  single  grid  domain  shown  in  Figure  3.3  with  uniform  spacing  in 
both  coordinate  directions.  Summing  the  above  equation  for  a typical  u control  volume  over 
all  the  u control  volumes  in  the  domain  yields 

nj  ni 

s = S X{(£e  - Ew)Ay  + (F"  - FAAx^ 

y=2 i=3 

nj  ni 

~ i _ ni  — Ew\i  - 3 )Ay  + 'y'{Fr^i  _ nj  — Fs\j  _ 2)Ax  = 0 
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Note  that  based  on  the  staggered  grid  convention  adopted  in  this  work,  this  procedure,  which 
is  taken  over  all  the  internal  unknown  u control  volumes,  results  in  a summation  over  the 
indices  i=3,  ni  and  j=2,  nj.  From  the  second  part  of  Equation  (3.3)  it  is  clear  that  S is  only 
a summation  of  the  boundary  fluxes,  since  the  interior  fluxes  for  neighboring  control 
volumes  cancel  each  other  out.  This  equation  represents  the  global  conservation  property  for 
any  control-volume  based  scheme  as  expressed  in  Equation  (3.2).  Notice,  however,  for  the 
staggered  grid,  the  boundary  in  question  is  not  the  physical  boundary  of  the  domain  as  is 
usually  the  case  when  dealing  with  a nonstaggered  grid,  when  the  unknowns  are  located  at 
the  cell  centers,  but  the  boundary  formed  by  the  boundary  sides  of  the  u control  volumes 
along  the  physical  boundary.  This  boundary  is  hereafter  referred  to  as  the  global 
conservation  boundary.  For  the  staggered  grid  system,  three  different  global  conservation 
boundaries  exist,  one  for  the  continuity  (pressure  correction)  equation,  one  for  the 
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(a) 


(c) 


Figure  3.4.  Global  conservation  boundaries  for  the  staggered  grid, 
(a)  Continuity;  (b)  u-momentum;  (c)  v-momentum. 


u-momentum  equation,  and  one  for  the  v-momentum  equation.  The  three  global 
conservation  boundaries  for  the  staggered  grid  covered  by  the  control  volumes  of  the  interior 
unknowns  are  shown  in  Figure  3.4. 

For  composite  grids,  the  global  conservation  property  also  requires  that  a summation 
of  the  u control  volume  equations  over  all  the  u control  volumes  results  in  only  a summation 
of  the  u momentum  fluxes  across  the  global  conservation  boundary  for  the  composite  grid. 
Note  that  for  such  a summation,  the  total  area  of  the  summation  should  be  equal  to  the  area 
formed  by  an  exclusive  summation  of  the  areas  enclosed  by  the  global  conservation 
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boundaries  of  each  of  the  blocks  forming  the  composite  grid.  By  an  exclusive  summation, 
we  mean  that  in  an  overlap  zone  created  by  the  intersection  of  the  global  conservation 
boundaries  of  two  blocks,  the  overlap  area  should  only  be  included  in  the  control  volume 
summation  of  one  of  the  two  blocks. 

For  example,  consider  the  composite  grid  shown  in  Figure  3.5(a),  constructed  from 
two  blocks  with  a horizontal  overlap  region.  The  global  conservation  boundary  for  the 
u-momentum  equation  is  highlighted  in  bold.  If  in  the  overlap  region,  we  sum  over  u control 
volumes  in  the  upper  block  only,  then  the  exclusive  overlapping  appears  as  shown  in  Figure 
3.5(b).  Numbers  have  been  given  to  the  various  boundary  segments  which  will  be  used  in 
the  ensuing  derivation  of  the  conditions  for  global  u-momentum  conservation.  Summing 
over  all  the  u control  volumes,  using  the  exclusive  overlapping  shown  in  Figure  3.5(b), 
yields 
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The  left-hand  side  of  this  expression  represents  the  summation  of  the  boundary  fluxes 
through  the  u-momentum  global  conservation  boundary  for  the  composite  grid.  The 
right-hand  side  of  the  expression  represents  the  difference  in  the  summation  of  the  fluxes 
along  the  internal  boundaries  of  the  upper  block  as  estimated  from  the  two  different  blocks 
in  the  composite  grid  system.  Now,  for  global  conservation  of  u-momentum,  the  summation 
of  the  u-momentum  flux  through  the  global  conservation  boundary  should  vanish. 


27 


Accordingly,  the  right-hand  side  of  Equation  (3.4)  should  vanish,  resulting  in  the  following 
necessary  condition  for  global  u-momentum  conservation 
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A more  restrictive  condition  can  be  imposed  which  requires  that  each  of  the  terms  in 
Equation  (3.5)  vanish  individually,  resulting  in 
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Figure  3.5.  Two-block  composite  grid  with  horizontal  overlap  region,  (a)  Composite 
grid  with  global  conservation  boundary  in  bold;  (b)  Exclusive 
overlapping  used  in  control  volume  summation. 
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These  relations  imply  that  the  u-momentum  flux  through  any  internal  boundary  must  be 
identical  when  estimated  from  the  blocks  on  either  side  of  the  internal  boundary.  It  should 
be  noted,  however,  that  these  are  only  sufficient  conditions  for  global  u-momentum 
conservation,  since  Equation  (3.5)  represents  the  only  requirement  for  global  u-momentum 
conservation.  Equation  (3.5)  states  that  for  global  conservation  of  u-momentum,  no  net 
u-momentum  can  be  created  at  the  internal  boundaries.  Similar  necessary  conditions  apply 
for  the  v-momentum  and  continuity  equations. 

3.3  Conservation  Strategy  for  Pressure-Based  Method 

Because  we  are  using  a pressure  correction  technique  for  solving  the  governing 
equations,  the  global  conservation  procedure  involves  a different  strategy  than  that  used  with 
density-based  compressible  flow  techniques  employing  nonstaggered  grids.  This  is  due  in 
part  to  the  nature  of  the  staggered  grid  system,  but  more  specifically  it  is  due  to  the  type  of 
boundary  condition  used  in  solving  the  pressure  correction  equation.  Two  types  of  boundary 
conditions  in  general  can  be  used  for  the  pressure  correction  equation.  If  the  pressure  is 
known  at  the  boundary,  then  the  value  of  the  pressure  correction  at  the  boundary  can  be  set 
to  zero.  If  the  pressure  is  not  known  at  the  boundary,  then  the  velocity  component  normal 
to  the  boundary  must  be  specified.  Since,  in  general,  the  boundary  pressure  values  are  not 
known,  especially  for  blocks  which  are  located  completely  interior  to  the  physical 
boundaries  of  the  domain,  we  have  adopted  the  use  of  the  normal  velocity  component 
boundary  condition  in  this  work  exclusively.  Normal  velocity  component  boundary 
conditions,  in  fact,  can  be  used  quite  generally  throughout,  even  at  boundaries  where  the 
normal  velocity  component  is  not  initially  known  (such  as  internal  boundaries  and  outflow 
boundaries),  but  which  evolves  to  a known  value  as  the  iterative  solution  process  converges. 
By  using  a normal  velocity  component  boundary  condition  for  the  pressure  correction 
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equation,  the  pressure  field  is  only  obtained  to  within  an  arbitrary  constant;  however,  this 
represents  no  problem  in  terms  of  the  solution  technique,  since  the  density  is  unaffected  by 
the  magnitude  of  the  pressure,  and  for  single  grid  solutions,  if  a pressure  value  is  known  at 
a certain  point  in  the  field,  then  the  entire  pressure  field  can  be  adjusted  accordingly,  after 
the  solution  has  been  obtained.  However,  for  composite  grids,  in  which  the  governing 
equations  are  solved  independently  within  each  block,  in  a block  to  block  iterative  fashion, 
each  block  in  the  system  will  generate  a pressure  field  independently  from  those  created  in 
neighboring  blocks.  It  is  this  characteristic  that  requires  a different  global  conservation 
strategy  than  that  used  in  density-based  compressible  flow  solvers  applied  to  composite 
grids. 

Solution  of  the  continuity  (pressure  correction)  equation  requires  the  specification 
of  the  mass  flux  into  each  of  the  pressure  correction  control  volumes  located  along  the 
boundaries  of  the  grid.  Similarly,  solution  of  the  u-  and  v-momentum  equations  requires  the 
specification  of  the  u-  and  v-momentum  fluxes  into  each  of  the  respective  u and  v control 
volumes  located  along  the  boundaries  of  the  grid.  Each  of  these  control  volume  fluxes  may 
involve  contributions  from  fluxes  due  to  specified  external  boundary  condition  segments, 
or  fluxes  entering  through  internal  boundaries  of  the  grid.  For  the  case  of  internal  boundaries, 
flux  contributions  must  be  calculated  in  such  a way  that  the  global  conservation  conditions 
previously  outlined  are  satisfied.  A discussion  of  the  exact  procedure  for  achieving  this  will 
be  undertaken  in  the  following  sections.  This  section  deals  with  issues  concerning  the  global 
conservation  strategy  that  has  been  developed  for  the  pressure  correction  algorithm  with  a 
staggered  grid  system. 

The  global  conservation  strategy  can  be  summarized  as  follows.  First,  explicit 
conservation  of  mass  across  the  horizontal  and  vertical  sides  of  the  global  conservation 
boundary  for  the  pressure  correction  equation  is  used  to  determine  the  mass  fluxes  into  each 
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of  the  pressure  correction  control  volumes  along  the  boundaries  of  the  grid  block  in  question. 
By  an  explicit  conservation  procedure,  we  mean  that  the  mass  flux  through  the  boundaries 
is  computed  directly  from  the  imposed  boundary  condition  segments  and  using  information 
from  the  neighboring  blocks  corresponding  to  the  internal  boundary  segments.  Once  the 
mass  fluxes  have  been  determined,  then  the  normal  component  of  the  velocity  profile  along 
the  horizontal  and  vertical  sides  of  the  grid  block  can  be  determined.  From  these  velocity 
values,  we  can  compute  the  contributions  to  the  u-momentum  fluxes  into  each  of  the  u 
control  volumes  along  the  vertical  boundaries  of  the  grid  block  which  do  not  involve 
pressure  (e.g.  the  term  guu  — //  ( du/dx  )),  as  well  as  the  contributions  to  the  v-momentum 
fluxes  into  each  of  the  v control  volumes  along  the  horizontal  boundaries  which  do  not 
involve  pressure  (e.g.  the  term  gw  — /u  (dv/ dy  )).  The  pressure  values  required  for 
calculating  the  rest  of  these  fluxes  are  already  specified  via  the  staggered  grid  arrangement. 
Since  the  u-momentum  fluxes  into  the  u control  volumes  along  the  horizontal  boundaries 
do  not  contain  a pressure  term  (i.e.  guv  - /i  ( du/dy  )),  they  can  be  obtained  via  explicit 
conservation,  as  was  done  for  the  mass  fluxes.  Similarly,  this  can  be  done  for  the 
v-momentum  fluxes  into  the  v control  volumes  along  the  vertical  boundaries.  In  this  way, 
all  the  required  boundary  fluxes  for  the  various  control  volumes  along  the  boundaries  of  the 
staggered  grid  block  in  question  are  obtained. 

For  a converged  solution  to  a composite  grid  problem,  since  the  pressure  fields 
within  each  of  the  blocks  of  the  system  have  developed  independently  from  the  others  and 
are  only  determined  to  within  a constant  within  each  grid  block,  the  pressure  fields  must  be 
corrected  in  some  manner  so  that  they  are  compatible  with  each  other.  In  order  to  determine 
the  compatibility  constant  between  two  blocks  sharing  an  internal  boundary,  a normal 
momentum  balance  is  performed  across  the  internal  boundary  interface,  and  the  pressure 
field  in  one  of  the  blocks  sharing  that  interface  is  adjusted  by  a constant  value  so  that  normal 


31 


momentum  flux  across  that  segment  is  identical  when  computed  from  either  of  the  two 
blocks.  For  the  general  case  in  which  the  grid  lines  from  neighboring  blocks  across  the 
interface  are  discontinuous,  this  post-processing  procedure  is  dependent  upon  the  choice  of 
the  segment  over  which  the  normal  momentum  flux  balance  is  performed;  however,  one  can 
compute  the  appropriate  compatibility  constants  for  ensuring  global  conservation  of  the 
normal  momentum  flux  by  performing  the  balance  described  above  along  an  appropriate 
length  of  the  internal  boundary.  This  issue  is  discussed  in  more  detail  in  a later  section. 

3.4  Explicit  Conservation  Procedure 

In  order  to  update  the  values  of  the  dependent  variables  in  a particular  block  of  a 
composite  grid,  the  fluxes  of  mass  and  momentum  through  the  boundaries  of  the  block  into 
the  various  boundary  control  volumes  must  be  specified.  According  to  the  global 
conservation  strategy  previously  outlined,  explicit  conservation  of  the  mass  and 
u-momentum  flux  through  the  horizontal  boundaries  and  the  mass  and  v-momentum  flux 
through  the  vertical  boundaries,  provides  all  of  the  information  required  for  calculating  the 
fluxes  into  all  of  the  boundary  control  volumes.  In  order  to  illustrate  the  explicit  conservation 
procedure  used  in  this  work,  the  explicit  conservation  of  mass  and  u-momentum  through 
horizontal  internal  boundaries  is  presented  in  detail  here.  The  same  concepts  apply  for 
explicitly  conserving  mass  and  v-momentum  through  vertical  internal  boundaries. 

Consider  a composite  grid  comprised  of  two  blocks  with  a horizontal  overlap  region. 
Figure  3.6  shows  a blowup  of  part  of  the  overlap  region,  where  the  upper  block  (whose 
dependent  variables  are  to  be  updated),  is  drawn  with  dotted  lines  and  the  neighboring  lower 
block  is  drawn  with  solid  lines.  The  upper  block  will  hereafter  be  denoted  as  block  one  and 
the  lower  block  as  block  two.  In  the  figure,  a pressure  correction  control  volume  along  the 
lower  boundary  of  block  one  has  been  shaded.  In  order  to  formulate  the  discrete  form  of  the 
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pressure  correction  equation  for  this  control  volume,  the  mass  flux  through  the  lower 
boundary  of  the  control  volume  (segment  AB)  must  be  obtained.  Since  segment  AB  lies 
completely  within  block  two,  all  the  information  required  for  calculating  the  mass  flux 
through  the  segment  is  obtained  from  this  block.  Segment  AB  is  comprised  of  a number  of 
sub-segments,  formed  by  the  intersection  of  the  vertical  grid  lines  of  block  two  with  the 
segment.  In  this  case,  segment  AB  is  comprised  of  five  complete  sub-segments  (forming  the 
center  portion  of  the  segment)  and  two  partial  segments  (on  the  left  and  right  ends  of  segment 
AB).  Each  of  these  sub-segments  contributes  a portion  to  the  total  mass  flux  through  AB. 

The  mass  flux  contribution  for  a typical  sub-segment  in  block  two  is  obtained  in  the 
following  manner.  The  constant  normal  velocity  component  along  the  segment  is  obtained 


Figure  3.6.  Blowup  of  horizontal  overlap  region  for  demonstrating  explicit 
conservation  of  mass.  Pressure  correction  control  volume 
is  shaded. 
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using  a linear  distance-weighted  interpolation  from  the  velocity  components  just  above  and 
below,  located  within  block  two.  For  example,  consider  the  second  sub-segment,  denoted 
as  ab,  of  segment  AB,  shown  in  Figure  3.6.  The  normal  velocity  component  along  the  entire 
sub-segment  is  taken  to  be  that  at  the  center  of  the  sub-segment  and  is  denoted  v2.  The  value 
of  v2  is  obtained  based  on  a linear  interpolation  within  block  two  from  the  staggered  grid 
values  located  directly  above  and  below,  which  are  also  shown  in  the  figure.  Once  the 
segment  velocity  has  been  found,  the  mass  flux  contribution  is  obtained  by  multiplying  this 
value  by  the  density  and  the  segment  length.  This  calculation  procedure  results  in  a piecewise 
constant  mass  flux  distribution  over  the  segment  AB.  With  the  total  mass  flux  from  block 
two,  denoted  as  mflux2,  and  assuming  that  the  north,  west,  and  east  pressure  correction 
neighbors  are  internal  unknowns  of  block  one,  the  discretized  form  of  the  pressure  correction 
equation  for  the  pressure  correction  control  volume  ABCD  becomes 

ap  P' p = aE  p' e aw  p'  w + aN  P' n + b (3.7a) 

aP  = aE  + aW  "k  aN  (3.7b) 

where  the  coefficients  a,-  are  as  given  before  (except  fora*  which  has  been  set  to  zero  since 
we  consider  the  mass  flux  obtained  to  be  a known  quantity)  and  the  source  term  b takes  the 
form 

b = [(gu*)w  - (Qu*)e](4y)i  + mflux2  - (gv*)n  (Ax)x  (3.7c) 

Once  the  mass  flux  into  the  pressure  correction  control  volume  has  been  calculated, 
the  velocity  component  at  the  control  volume  boundary  for  block  one,  denoted  vb  in  the 
figure,  is  given  by  the  following: 

mflux2 

Vb  ~ 


(3.8) 
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Thus,  as  previously  stated,  conservation  of  mass  along  the  horizontal  boundaries  provides 
both  the  boundary  condition  for  each  of  the  pressure  correction  control  volumes  along  the 
boundary,  as  well  as  the  normal  velocity  component  profile  along  the  boundary.  Once  this 
profile  is  determined,  then  all  the  information  required  for  the  calculation  of  the 
contributions  to  the  v-momentum  fluxes  not  involving  pressure  (i.e.  qw  — /x  (dv/dy))is 
available. 

The  procedure  for  explicitly  conserving  the  u-momentum  fluxes  into  each  of  the  u 
control  volumes  along  horizontal  boundaries  is  similar  to  that  for  explicitly  conserving  mass. 
Consider  again  a composite  grid  comprised  of  two  blocks  with  a horizontal  overlap  region. 
Figure  3.7  shows  a blowup  of  part  of  the  overlap  region,  where  the  upper  block  (whose 
dependent  variables  are  to  be  updated)  is  again  denoted  as  block  one  and  the  lower 
neighboring  block  as  block  two.  In  the  figure,  a u control  volume,  labelled  ABCD,  along  the 
lower  boundary  of  the  block  has  been  shaded.  In  order  to  formulate  the  discrete  form  of  the 
u-momentum  equation  for  this  control  volume,  the  u-momentum  flux  through  the  lower 
boundary  of  the  control  volume  (segment  AB)  must  be  obtained.  Again,  segment  AB  is 
comprised  of  a number  of  sub-segments;  however,  in  this  case  the  sub-segments  are  defined 
by  the  u control  volumes  in  block  two.  Consider  the  sub-segment  ab  shown  in  the  figure.  The 
u-momentum  flux,  including  the  convective  and  diffusive  components,  through  this 
sub-segment,  denoted  ufluxab  is  calculated  as  follows: 

ufluxab  = [D(l  - 0.5IF/DI)  + |[F,  0.01  ] uL  - (3.9a) 

[D(l  - 0.5IF/DI)  + [[  - F,  0.0J  ] uv 

where  the  terms  Uu  and  uL  are  respectively  the  u velocity  components  located  on  the 
staggered  grid  of  block  two  just  above  and  below  the  segment.  In  this  expression,  the 
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second-order  central  difference  scheme  has  been  used  for  both  the  convection  and  diffusion 
terms  for  illustration  purposes.  The  terms  F and  D represent  to  the  flux  trough  the  segment 
from  convection  and  diffusion  and  are  given  by 

D = ^7  F = CV.0X)  (3.9b) 

where  Zlx  represents  the  length  of  the  sub-segment  and  (Ay)2  is  the  vertical  distance  between 
uL  and  u jj  which  for  this  case  is  equal  to  the  constant  vertical  grid  spacing  in  block  two.  The 
quantity  vs  is  the  velocity  component  calculated  via  bi-linear  interpolation  from  the 
neighboring  v component  values  of  block  two  at  the  center  of  the  segment.  With  the  total 
flux  contribution  from  block  two  denoted  as  uflux2,  and  assuming  that  the  north, east,  and 
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west  neighbors  of  the  u boundary  control  volume  in  block  one  are  internal  unknowns,  the 
discretized  form  of  the  u-momentum  equation  for  the  control  volume  becomes 


In  the  original  formulation  for  the  u control  volume  shown,  the  term  (Fn  - Fs  + Fe  - Fw ) 
in  Equation  (2.3e),  is  identically  zero  for  a continuity  satisfying  velocity  field  and  can  be 
dropped.  The  term  ( Fn  + Fw  - Fe)uP  remains  in  Equation  (3.10a)  because  AP  is  now 
computed  without  the  contribution  from  As.  Due  to  the  interface  treatment,  As  and  its 
corresponding  component  in  AP  are  now  explicitly  given  in  the  form  of  uflux2.  The  term 
(Fn  + Fw  - Fe)  can  be  replaced  by  Fs  and  computed  using  the  normal  velocity  component 
profile  at  the  boundary  obtained  from  explicitly  conserving  mass. 

In  the  examples  above,  illustrating  the  explicit  procedure  for  conserving  mass  and 
u-momentum  across  horizontal  boundaries,  all  the  required  information  was  obtained  from 
a single  neighboring  block.  It  should  be  noted  here,  that  in  a general  composite  grid,  for  any 
control  volume  located  along  a boundary,  contributions  to  the  total  flux  into  the  control 
volume  may  be  required  from  any  number  of  neighboring  blocks  and/or  externally  imposed 
boundary  conditions.  Since  the  path  of  information  flow  into  each  block  is  entirely  specified 
via  boundary  segments,  we  can  determine  which  neighboring  block  or  boundary  condition 
is  to  provide  the  required  information  across  a particular  portion  of  a control  volume 
boundary.  Once  this  is  known,  then  if  a neighboring  block  is  required  to  provide  the 
information,  the  flux  through  that  portion  can  be  calculated  as  described  above,  and  if  a 


Ap  UP  — Ag  Mg  + A]y  Uyy  + Ajy  Upj  — (F n + Fe  — Fw)uP  + 
(Pw  ~ PE)(Ay)x  + uflux2 


(3.10a) 


Ap  — AE  + Ayy  + AN 


(3.10b) 


37 


boundary  condition  is  required,  then  the  flux  across  that  portion  can  be  calculated 
accordingly. 


3.5  Conditions  for  Local  and  Global  Conservation 

In  this  work,  as  previously  mentioned,  the  only  restriction  placed  on  the  overlap  of 
grid  blocks  is  that  internal  block  boundaries  must  fall  within  a constant  column  or  row  index 
of  any  neighboring  blocks.  For  such  configurations,  the  question  arises  as  to  whether  both 
local  and  global  conservation  can  be  maintained  simultaneously.  Alternatively  posed,  does 
the  local  conservation  procedure  detailed  in  the  previous  section  result  in  a globally 
conservative  formulation  for  general  configurations?  In  accordance  with  Patankar’s 
requirement  that  acceptable  formulations  apply  equally  well  to  fine  and  coarse  grid 
discretizations,  a similar  requirement  is  placed  here  on  the  conservative  interface  treatment, 
namely  that  we  be  able  to  drive  global  mass  and  momentum  residuals  to  machine  accuracy 
for  arbitrary  configurations  satisfying  the  single  constraint  mentioned  above. 

To  begin  the  discussion,  consider  the  following  two-block  cavity  flow  configuration 
shown  in  Figure  3.8.  For  the  purpose  of  clarity,  the  two  blocks  are  shown  separated;  however, 
physically  they  overlap.  The  dotted  vertical  line  shown  within  the  left  block  (block  one), 
represents  the  physical  location  where  the  left  boundary  of  the  right  block  (block  two)  falls 
within  block  one.  Suppose  that  the  mass  fluxes  into  the  control  volumes  along  the  left 
boundary  of  block  two  are  determined  from  block  one  using  the  explicit  conservation 
procedure  previously  detailed.  In  order  to  satisfy  global  mass  conservation  for  block  two, 
the  mass  fluxes  must  sum  to  zero,  since  no  mass  can  enter  through  the  upper,  lower,  and  right 
solid  boundaries.  From  the  explicit  conservation  of  mass  into  block  two,  the  total  mass  flux 
through  line  C,  denoted  as  Mc,  can  be  written  in  terms  of  the  mass  fluxes  through  lines  A 


and  B as 
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where  xA  and  xB  represent  the  horizontal  coordinates  of  lines  A and  B respectively  and  Jtthe 
horizontal  coordinate  of  the  dotted  line.  Since  we  are  using  a control  volume  formulation, 
both  Ma  and  MB  are  identically  zero  (to  machine  precision),  and  thus,  Mc  = 0 holds. 
Therefore,  for  this  case,  the  explicit  local  conservation  procedure  also  results  in  a globally 
conservative  formulation  for  mass.  One  important  thing  to  note  here  is  that  due  to  the 
piecewise  constant  local  conservation  procedure,  the  total  fluxes  MA  and  MB  are  consistent 
with  the  internal  mass  summations  along  the  lines  A and  B and  are  thus  by  definition  zero. 
If  a piecewise  linear  local  conservation  (in  the  vertical  direction)  is  used,  this  is  not  the  case 
and  in  general  a global  mass  correction  must  be  distributed  along  line  C to  enforce  global 
conservation,  which  destroys  any  concept  of  local  conservation.  Another  interesting  thing 
to  note  is  that  as  far  as  v-momentum  is  concerned,  for  this  case,  global  conservation  is  always 
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maintained,  regardless  of  the  explicit  conservation  procedure  used  along  line  C,  due  to  the 
presence  of  the  solid  boundaries,  which  will  adjust  accordingly  to  produce  a net 
v-momentum  flux  of  zero  from  block  two.  In  this  regard,  solid  walls,  in  general,  provide  an 
additional  degree  of  freedom  for  ensuring  global  conservation  of  momentum;  however,  they 
allow  no  such  flexibility  as  far  as  mass  is  concerned.  Outflow  boundaries,  on  the  other  hand, 
can  provide  relief  for  both  mass  and  momentum.  For  grid  blocks  which  are  located  entirely 
internal  to  a domain;  however,  a global  imbalance  of  momentum  can  occur  as  with  mass 
above,  since  no  additional  degree  of  freedom  is  provided  in  this  case  via  either  a solid  wall 
or  outflow  condition. 

Now  consider  the  three-block  cavity  flow  configuration  shown  in  Figure  3.9.  Here 
we  are  concerned  with  the  total  mass  flux  entering  block  three.  For  this  case,  the  overlapping 
priorities  are  set  with  block  one  having  the  highest  priority,  followed  by  block  two  and  then 
block  three.  With  this  arrangement,  the  internal  boundaries  of  block  three  are  divided  into 
three  segments,  labelled  A,  B,  and  C at  the  center  of  the  segments,  as  shown.  The  total  mass 
flux  entering  through  the  internal  boundaries  of  block  three  is  given  as 

mflux3  = Ma  + Mb  + Mc  (3.12) 

For  global  conservation,  mflux3  must  be  identically  zero,  since  no  mass  can  enter  block  three 
through  the  lower  or  right  solid  walls.  In  this  case,  however,  unlike  the  previous  one,  no 
control  volume-based  argument  can  be  made  to  prove  that  this  summation  is  zero,  and  in 
general,  the  summation  is  not  zero.  Therefore,  a global  mass  imbalance  for  block  three 
results,  preventing  the  continuity  (pressure  correction)  equation  from  driving  its  residual  to 
machine  zero.  The  level  of  mass  imbalance  for  this  scenario  is  also  highly  dependent  on  both 
the  discretization  as  well  as  the  solution  gradients  within  the  domain,  since  both  of  these 
factors  influence  the  accuracy  of  the  local  interpolation.  For  a sufficiently  fine  discretization, 
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global  mass  conservation  may  be  closely  satisfied,  but  only  in  the  limit  of  infinite  grid 
resolution  is  it  guaranteed  to  be  exactly  satisfied.  Again  though,  no  global  conservation 
problem  exists  for  the  u-momentum  and  v-momentum  equations  due  the  presence  of  the 
solid  walls. 

In  general,  in  order  to  achieve  a globally  conservative  formulation  via  the  piecewise 
constant  explicit  local  conservation,  all  locally  conserved  fluxes  must  be  based  strictly  upon 
the  interior  fluxes  of  neighboring  blocks  as  they  are  represented  within  the  control  volume 
formulation,  without  resorting  to  interpolation  among  the  fluxes.  For  the  global  conservation 
of  mass,  this  implies  that  the  internal  boundaries  of  every  block  in  the  composite  grid  must 
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Figure  3.9.  Three-block  cavity  flow  configuration.  Explicit  conservation  procedure 
does  not  guarantee  global  conservation  of  mass  for  block  three. 

Internal  boundaries  are  highlighted  in  bold. 
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fall  identically  upon  the  grid  lines  of  the  neighboring  blocks.  Since  mass  and  tangential 
momentum  are  the  only  quantities  which  are  explicitly  conserved,  and  since  they  share  the 
same  physical  internal  boundary  on  every  side  of  a block,  this,  in  fact,  is  the  only  condition 
required  to  guarantee  global  conservation  of  mass  and  both  momentum  components  for 
general  configurations.  This  overlap  condition  for  ensuring  global  conservation  is  similar 
to  that  imposed  by  Rai  (1986a)  in  his  development  of  a conservative  transfer  scheme  for 
patched  grid,  compressible  flows,  where  artificial  grid  lines  are  extended  into  the 
neighboring  blocks  until  they  intersect  the  nearest  interior  control  volume  face  at  the 
boundary,  which  in  effect  results  in  a direct  usage  of  the  neighboring  block  fluxes  (with  the 
piecewise  constant  local  conservation)  without  resorting  to  interpolation. 

3.6  Numerical  Validation:  Lid-Driven  Cavity  Flow 

To  assess  the  functionality  of  the  organizational  procedure  and  to  provide  a 
numerical  validation  of  the  interfacing  scheme,  the  flow  in  a square  lid-driven  cavity  is 
computed  using  the  composite  grid  approach.  This  flow  is  commonly  accepted  as  a standard 
of  comparison  for  the  solution  of  the  incompressible  Navier-Stokes  equations.  For 
comparison  purposes  here,  benchmark  solutions  obtained  by  Ghia  et  al.  (1982)  are  used.  The 
flows  for  two  Reynolds  numbers,  Re=400  and  Re=1000,  are  computed.  These  values  are 
based  on  the  sliding  lid  velocity  and  width  of  the  domain.  For  the  computations  shown  here, 
the  width  of  the  domain  is  taken  to  be  unity  and  the  infinite  lid  is  prescribed  to  be  sliding  to 
the  left  with  a velocity  of  unity.  The  block  partitioning  and  grid  used  for  both  computations 
are  shown  in  Figure  3.10.  The  grid  is  composed  of  nine  individual  grid  blocks,  each  with 
either  a constant  grid  spacing  of  1/80  or  1/160.  Finer  grid  blocks  have  been  placed  near 
regions  of  anticipated  secondary  recirculation  and  in  regions  where  the  flowfield  gradients 
are  anticipated  to  be  the  sharpest.  It  should  be  noted  that  the  grid  used  here  was  not 
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Figure  3.10.  Lid-driven  cavity  flow,  (a)  Block  partitioning;  (b)  Grid. 
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constructed  as  an  optimal  grid  for  the  particular  flows  being  computed,  but  instead  was 
created  mainly  to  demonstrate  the  capability  of  the  organizational  strategy  while  at  the  same 
time  verifying  the  ability  of  the  interface  scheme  to  correctly  handle  grid  overlaps  in  which 
discontinuous  grid  lines  occur  across  internal  interfaces. 

Streamlines  for  the  two  computed  cases  are  shown  in  Figure  3.11.  For  both  flows, 
second-order  central  differences  were  used  for  the  convection  and  diffusion  terms. 
Physically,  both  flows  consist  of  a primary  recirculation  region  centered  roughly  just  above 
the  geometric  center  of  the  cavity,  along  with  secondary  recirculation  zones  in  each  of  the 
lower  comers.  As  the  Reynolds  number  is  increased  from  400  to  1000,  the  center  of  the 
primary  recirculation  region  shifts  downward,  towards  the  geometric  center  of  the  cavity. 
In  addition,  the  secondary  recirculation  zones,  most  notably  the  one  in  the  lower  right  comer, 
increase  in  size.  Both  of  these  qualitative  assessments  were  also  found  by  Ghia  et  al.  (1982). 
Due  to  the  conservative  treatment  of  the  mass  flux  at  the  internal  boundaries,  no 
inconsistency  in  the  streamfunction  occurs  across  the  internal  boundaries  and  both  the 
primary  and  secondary  recirculation  zones  pass  smoothly  through  the  internal  boundaries 
without  distortion. 

In  terms  of  a quantitative  solution  comparison,  plots  of  the  u-velocity  component 
distribution  along  the  vertical  line  through  the  geometric  center  of  the  cavity,  and  the 
v- velocity  component  distribution  along  the  horizontal  line  through  the  geometric  center  are 
shown  in  Figure  3.12  for  both  flows.  The  benchmark  solutions  plotted  in  the  figure  were 
obtained  using  a 129x129  uniform  grid  along  with  a second-order  upwind  scheme  for  the 
convection  terms.  Based  on  the  total  number  of  nodes  used,  the  composite  grid  corresponds 
approximately  to  a 125x125  uniform  grid  resolution.  For  both  Reynolds  numbers,  our 
solutions  agree  very  well  with  the  benchmark. 
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Figure  3.11.  Streamlines  for  lid-driven  cavity  flow,  a)  Re=400;  b)  Re=1000. 
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3.7  Computation  in  Multiply-Connected  Flow  Domains 

With  the  present  interface  treatment,  since  the  mass  flux  provides  both  the  boundary 
condition  for  the  pressure  correction  equation  as  well  as  the  velocity  component  for  the 
normal  momentum  flux  at  the  boundary,  if  the  mass  flux  into  the  boundaries  of  each  block 
of  the  system  are  uniquely  obtained  during  the  course  of  the  calculation  procedure,  then  the 
resulting  pressure  fields  within  the  different  blocks  will  be  compatible  with  each  other,  up 
to  an  additive  constant.  Such  is  the  case  for  all  simply-connected  domains,  and  thus,  no 
exchange  of  normal  momentum  flux  across  the  internal  boundaries  is  required  during  the 
solution  process.  For  general  multiply-connected  regions;  however,  an  ambiguity  in  the 
mass  flux  distribution  among  the  blocks  can  occur  unless  the  normal  momentum  flux 
conservation  constraint  is  explicitly  maintained.  To  examine  this  issue  in  more  detail, 
consider  the  multiply-connected  flow  domain  shown  in  Figure  3.13,  consisting  of  a channel 
with  internal  baffles.  The  region  has  been  decomposed  into  seven  individual  computational 
zones,  with  the  overlap  regions  shown  by  dotted  lines. 

Consider  the  mass  flux  entering  and  exiting  block  two.  The  mass  flux  entering  block 
two  from  block  one  is  uniquely  specified  since  the  flux  imposed  at  the  inlet  of  block  one 
completely  enters  block  two.  During  the  course  of  the  iterative  solution  process,  this  total 
flux  is  distributed  in  some  manner  to  provide  the  inlet  conditions  for  blocks  three,  four,  and 
five.  For  a converged  solution  with  the  appropriate  mass  distribution  among  blocks  three, 
four,  and  five,  the  pressure  fields  within  all  the  blocks,  upon  adjustment  by  the  compatibility 
constants  will  be  consistent,  and  normal  momentum  flux  will  be  conserved  across  all  the 
internal  boundaries.  Since,  however,  no  pressure  information,  and  thus,  no  normal 
momentum  flux  is  exchanged  between  the  blocks  during  the  solution  process,  any  arbitrary 
mass  distribution  to  blocks  three,  four  and  five  which  satisfies  the  mass  conservation 
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constraint  is  acceptable,  and  thus,  the  possibility  of  multiple  solutions  exists.  In  fact,  the  mass 
flux  distribution  among  blocks  three,  four,  and  five  can  vary  arbitrarily  with  the  overlap 
thickness,  initial  guess  and  relaxation  parameters,  as  well  as  other  factors. 

In  order  to  test  the  compatibility  of  the  pressure  fields  of  the  various  blocks,  a 
pressure  cycling  sequence  must  be  established.  By  a pressure  cycling  sequence,  we  mean  the 
sequence  of  normal  momentum  flux  comparisons  between  blocks  used  to  determine  the 
compatibility  constants  for  each  of  the  blocks  in  the  composite  grid.  For  the  configuration 
shown  in  Figure  3. 13,  a reasonable  sequence  might  be  as  follows.  First  a normal  momentum 
flux  balance  is  performed  along  the  interface  separating  blocks  one  and  two  by  adjusting  the 
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Figure  3.13.  Channel  with  internal  baffles. 


48 


pressures  within  block  two  by  the  appropriate  constant.  Once  the  pressures  within  block  two 
have  been  adjusted,  then  similar  balances  can  be  performed  to  make  blocks  three,  four  and 
five  compatible  with  block  two.  For  block  six,  the  normal  momentum  balance  can  be 
performed  with  either  block  three,  four,  or  five;  however,  here  is  where  the  issue  of 
compatibility  arises.  Since  the  pressures  is  block  six  are  to  be  adjusted  only  once,  if  the 
normal  momentum  balance  is  performed  with  respect  to  one  of  the  blocks  (say  block  three), 
then  for  a unique  solution  in  which  the  mass  fluxes  have  been  distributed  correctly  among 
the  passageways,  the  pressure  fields  at  the  interfaces  of  the  other  two  blocks  (blocks  four  and 
five)  and  block  six  will  be  consistent.  In  other  words,  for  a unique  solution,  which  honors 
both  mass  and  normal  momentum  conservation,  the  compatibility  constants  between  blocks 
three  and  six,  blocks  four  and  six,  and  blocks  five  and  six,  are  identical,  and  thus,  the  pressure 
field  within  block  six  can  be  corrected  with  respect  to  any  one  of  the  three  blocks.  After  the 
pressure  field  within  block  six  has  been  corrected,  then  a normal  momentum  balance  at  the 
interface  with  block  seven  can  be  conducted  to  adjust  the  pressure  field  in  block  seven,  which 
completes  the  pressure  cycling  process. 

Without  explicitly  enforcing  the  normal  momentum  flux  conservation  constraint 
across  the  internal  boundaries,  upon  adjustment  of  the  pressure  field  in  block  six  with  respect 
to  block  three,  discontinuities  will  appear  at  the  interfaces  between  blocks  four  and  six  and 
blocks  five  and  six.  In  order  to  ensure  that  these  unphysical  jumps  do  not  occur,  the  following 
procedure  is  employed.  At  the  end  of  every  outer  iteration  (defined  as  a sweep  through  the 
system  of  equations  for  every  block),  the  pressure  cycling  procedure  detailed  above  is 
carried  out.  This  results  in  a smooth  pressure  field  across  all  internal  boundaries  except  those 
between  blocks  four  and  six  and  blocks  five  and  six.  After  the  pressure  cycling  is  completed, 
additional  normal  momentum  balances  are  performed  between  blocks  four  and  six  and 
blocks  five  and  six;  however,  only  the  pressure  values  in  the  first  row  of  cells  of  block  six 
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in  the  overlap  regions  with  blocks  four  and  five  are  adjusted.  This  results  in  pressure 
discontinuities  being  introduced  into  block  six  along  those  lines;  however,  it  is  these  pressure 
discontinuities  which  drive  the  overall  mass  flux  partitioning  between  the  three  passageways 
to  the  correct  values  which  satisfy  both  mass  conservation  and  normal  momentum 
conservation  across  all  internal  boundaries.  As  the  solution  progresses,  these  discontinuities 
become  diminished  and  ultimately  vanish  at  convergence. 

In  order  to  demonstrate  that  the  above  procedure  results  in  a physically  correct 
solution,  the  flow  in  the  configuration  shown  previously  in  Figure  3.13  is  computed  for  a 
Reynolds  number  of  of  800  based  on  the  entrance  velocity  and  the  configuration  width.  For 
simplicity,  a grid  with  continuous  grid  lines  across  the  internal  interfaces  has  been 
constructed,  as  shown  in  Figure  3.14.  The  resolutions  for  each  of  the  various  blocks  are  also 
shown  in  the  figure.  These  grid  resolutions  result  in  horizontal  overlaps  two  cells  thick 
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Figure  3.14.  Grid  used  to  compute  flow  shown  in  Figure  3.13. 
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between  blocks  one  and  two  and  blocks  six  and  seven,  and  five  cells  thick  between  blocks 
three,  four,  and  five  and  both  blocks  two  and  six.  Central  differences  have  been  used  here 
for  both  the  convection  and  diffusion  terms.  Streamfunction  and  pressure  contours  are 
presented  in  Figures  3.15(a)  and  3.15(b)  respectively.  Two  large  recirculation  regions  are 
found  to  exist  on  either  side  of  the  entrance  to  the  main  portion  of  the  channel.  A large 
recirculation  region  is  also  present  in  the  center  of  the  middle  passage  of  the  channel,  created 
due  to  the  large  flow  deflection  angle  caused  by  the  impingement  of  the  main  flow  on  the 
baffle  located  near  the  entrance.  Some  of  the  main  flow  which  initially  passed  through  the 
left-most  passage  of  the  channel  is  drawn  deep  into  the  middle  passage,  encompassing  the 
recirculation  region  there,  before  being  drawn  out  by  the  main  flow  through  the  middle 
passage  to  the  channel  exit.  From  the  pressure  contours,  it  is  clear  that  the  current  procedure 
does  produce  a physically  realistic  solution  with  no  artificial  pressure  jumps  across  the 
internal  interfaces. 

3.8  Implementation  of  Higher-Order  Convection  Schemes 

The  preceding  development  of  the  internal  boundary  strategy  has  concentrated  on  the 
use  of  standard  five-point  stencil  schemes  involving  two-point  approximations  for  the 
convection  terms  at  each  of  the  control  volume  faces.  Such  schemes  include  the  central 
difference,  hybrid,  and  first-order  upwind  schemes.  For  these  schemes,  conservation  of  mass 
at  the  internal  boundaries  provides  all  the  required  information  to  obtain  the  normal 
momentum  fluxes  since  only  the  normal  velocity  component  at  the  boundary  is  required, 
along  with  the  nearest  internal  neighbor,  to  compute  the  convection  term  at  the  boundary  for 
the  normal  momentum  control  volume.  In  order  to  obtain  better  accuracy  for  incompressible 
flows  and  compressible  flows  with  smooth  (shock-free)  solutions,  higher-order  convection 
schemes  such  as  the  second-order  upwind  scheme  and  the  QUICK  (Quadratic  Upstream 
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Figure  3.15.  Flow  through  channel  with  internal  baffles  (Re  = 800). 

(a)  Streamfunction  contours;  (b)  Pressure  contours. 
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Figure  3.15  — continued 


53 


Interpolation  for  Convection  Kinematics)  scheme  (Leonard  1979)  are  now  commonly 
employed.  Several  different  implementations  of  each  of  these  schemes  are  possible; 
however,  here  we  are  concerned  only  with  the  conservative  implementations  of  each  scheme, 
as  detailed  in  Thakur  (1993). 

To  illustrate  the  implementation  of  the  internal  boundary  scheme  for  such 
higher-order  convection  schemes,  the  second-order  upwind  scheme  will  be  used  as  the  basis 
for  discussion.  The  implementation  of  other  higher-order  schemes  follows  directly.  To 
review  the  formulation  of  a typical  control  volume  interface  flux  for  the  second-order 
upwind  scheme,  consider  the  interior  u control  volume  shown  in  Figure  3.16,  along  with  the 
two  upwind  neighbors  uw  and  mw.  For  this  case,  it  is  assumed  that  the  flow  comes  from 
the  left  direction. 


The  u-momentum  flux  at  the  west  control  volume  face  is  written  as 

(guu)wAy  + p^Ay  - (p^ Ay  (3.13) 

For  the  total  flux  shown  in  Equation  (3.13),  the  pressure  term  is  obtained  directly  due  to  the 
staggered  grid  and  the  diffusion  term  is  computed  via  a standard  two-point  central 
difference.  With  the  second-order  upwind  scheme,  the  convection  term  is  handled  in  the 
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following  manner.  The  convection  term  (gu^wAy  is  interpreted  as  {gu^wAy  uw  , where 
{gu)wAy  is  the  mass  flux  at  the  west  control  volume  face  and  is  obtained  from  the  velocity 
components  on  either  side  of  the  face.  For  a uniform  grid  spacing,  this  results  in 

(pwMy  = | (uw  + Up) Ay  (3.14) 

To  obtain  the  second-order  formulation,  the  dependent  variable  at  the  face,  uw , is  obtained 
via  a linear  extrapolation  from  the  upwind  neighbors.  For  a uniform  grid  spacing,  this  gives 

3 1 

uw  = 2uw  ~ 2 Mww  (3.15) 

and  thus  the  entire  convection  term  is  expressed  as 

Q 3 

2^uw  + uP){-^uw  — 2 uww)  (3.16) 

It  should  be  noted  here,  that  in  order  to  obtain  a conservative  form  for  the  convection  term, 
the  upwind  direction  must  be  evaluated  based  upon  the  sign  of  the  u-velocity  component 
value  at  the  control  volume  face,  which  is  indicated  by  the  mass  flux  term. 

With  the  above  discussion  in  mind,  the  internal  boundary  scheme  can  be  extended 
for  the  second-order  upwind  scheme  as  follows.  For  the  explicit  conservation  of  mass  at 
internal  boundaries  for  the  pressure  correction  control  volumes,  the  same  procedure  as 
detailed  in  Section  3.4  applies,  since  only  the  convection  terms  of  the  momentum  equations 
are  being  treated  differently  now.  Since  the  tangential  momentum  flux  through  the  internal 
boundaries  is  also  explicitly  conserved,  these  fluxes  can  be  computed  from  the  neighboring 
block  velocity  values  in  a manner  consistent  with  the  internal  representation  of  these  fluxes, 
as  described  above.  Therefore,  for  mass  and  tangential  momentum,  no  fundamental 
additions  need  to  be  added  to  the  current  interface  scheme.  As  far  as  the  normal  momentum 
component  is  concerned,  some  means  must  be  established  to  provide  an  additional  normal 
velocity  component  value  at  the  internal  boundaries,  so  that  the  convection  term  can  be 
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appropriately  handled.  Consider  a u control  volume  located  along  an  internal  boundary,  as 
shown  in  Figure  3.17.  If  the  flow  comes  from  the  left  direction,  the  velocity  components 
along  lines  A and  B,  as  well  as  the  nodal  value  itself  are  required  to  compute  the  convection 
flux  term.  Line  B corresponds  to  the  physical  grid  boundary,  and  thus,  the  velocity 
component  there  is  already  known  from  conservation  of  mass  for  the  pressure  correction 
control  volumes.  Line  A represents  an  imaginary  extension  from  the  block  under 
consideration  and  is  used  to  obtain  the  additional  required  component.  The  velocities  along 
line  A can  be  obtained  in  a manner  identical  to  that  used  for  line  B from  the  neighboring  block 
information.  To  ensure  that  this  procedure  for  computing  the  convection  term  degenerates 
appropriately  to  the  case  where  two  blocks  of  equal  grid  spacing  overlap  with  continuous 
grid  lines  across  the  interface,  the  spacing  between  lines  A and  B should  be  the  same  as  the 
internal  grid  spacing,  as  shown  in  the  figure. 

3.9  Concluding  Remarks 

In  this  chapter,  a detailed  discussion  of  a pressure-based  composite  grid  method  has 
been  presented.  An  organizational  scheme  has  been  developed  based  on  the  idea  of  boundary 
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Figure  3.17.  U control  volume  at  grid  boundary,  with  velocity  components 
required  for  second-order  upwind  scheme. 
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segments  which  allows  generic  configurations  composed  of  any  number  of  arbitrarily 
overlapping  blocks  to  be  handled  in  a straightforward  manner.  Such  a scheme  is  essential  for 
developing  a conservative  internal  boundary  treatment  for  arbitrary  flow  domains.  From  an 
examination  of  the  global  conservation  conditions  for  the  staggered  grid,  three  global 
conservation  boundaries  have  been  identified,  one  associated  with  each  of  the  continuity, 
u-momentum  and  v-momentum  equations.  The  existence  of  three  global  conservation 
boundaries  is  a defining  characteristic  of  the  staggered  grid  and  represents  a fundamental 
difference  from  the  so-called  collocated  grids  which  maintain  one  global  conservation 
boundary  for  all  unknowns. 

An  internal  boundary  scheme,  based  on  the  local  conservation  of  mass  and  tangential 
momentum  flux  has  been  proposed,  which  remains  true  to  the  fundamental  nature  of  the 
staggered  grid  arrangement,  in  which  no  artificial  pressure  boundary  conditions  are 
required.  For  simply-connected  domains,  global  conservation  of  the  normal  momentum  flux 
across  internal  boundaries  provides  the  means  for  eliminating  the  arbitrary  constants  of  the 
pressure  field  which  develop  due  to  the  decoupling  nature  of  the  normal  velocity  component 
boundary  condition  used  for  the  pressure  correction  equation.  For  multiply-connected  flow 
domains,  global  conservation  of  normal  momentum  flux,  within  the  context  of  a 
pressure-cycling  strategy,  has  been  shown  to  play  a more  fundamental  role  in  determining 
unique  solutions  in  which  all  flux  quantities  are  conserved  across  internal  boundaries. 
Finally,  extension  of  the  current  methodology  to  incorporate  higher-order  convection 
schemes  requires  no  fundamental  modifications  to  the  existing  interface  strategy. 


CHAPTER  4 

COMPOSITE  MULTIGRID 
4.1  Background 

In  order  to  substantially  increase  the  computational  efficiency  of  the  current 
capability,  a multigrid  procedure,  previously  developed  by  Shyy  and  Sun  (1993)  for 
single-block  staggered  grid  domains,  has  been  extended  to  the  present  composite  grid 
algorithm.  The  motivation  behind  the  development  of  the  multigrid  (or  multilevel)  method 
comes  from  the  realization  that  for  typical  relaxation  techniques,  fast  convergence  can  only 
be  achieved  in  the  initial  iterations.  By  considering  a Fourier  analysis  of  the  error  reduction 
process  for  typical  relaxation  methods,  Brandt  (1977)  showed  that  these  techniques  are  only 
efficient  at  eliminating  error  components  whose  wavelengths  are  comparable  to  the  finite 
difference  mesh  spacing.  These  components  are  quickly  eliminated  by  the  iterative  process, 
while  longer  wavelength  error  components  are  eliminated  much  more  slowly.  Thus, 
relaxation  techniques  have  also  come  to  be  commonly  known  as  ’’smoothing”  techniques, 
since  only  the  long  wavelength  components  remain  after  a few  iterations. 

The  fact  that  the  rate  of  convergence  decreases  sharply  after  the  first  few  iterations 
represents  a fundamental  problem  in  terms  of  efficiently  solving  systems  of  equations  via 
relaxation  techniques.  In  many  instances,  where  a large  number  of  grid  points  are  employed 
in  the  discretization,  relaxation  techniques  offer  the  only  viable  alternative  for  solving  the 
resulting  system  of  linear  equations,  as  the  direct  inversion  of  the  coefficient  matrix  becomes 
intractable  due  to  computer  memory  limitations.  However,  because  the  rate  of  convergence 
is  tied  directly  to  the  spectral  radius  of  the  coefficient  matrix,  which  typically  approaches 
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unity  as  the  mesh  spacing  is  reduced,  convergence  problems  become  even  more  severe  for 
cases  employing  very  fine  grids  such  as  those  used  in  obtaining  accurate  solutions  for 
complex  fluid  flows.  Since  the  CPU  time  also  increases  in  direct  relation  to  the  number  of 
grid  points,  a prohibitive  amount  of  CPU  time  may  be  required  for  solving  large  problems 
with  the  use  of  standard  relaxation  methods  alone. 

The  multigrid  method  recognizes  that  an  error  component  is  most  efficiently 
eliminated  on  a grid  whose  mesh  size  is  comparable  to  its  wavelength,  and  thus  a sequence 
of  grids  of  different  mesh  spacing  is  used  to  eliminate  different  components  of  the  overall 
error  in  the  most  efficient  manner.  This  sequence  of  grids  will  hereafter  be  denoted  by  Gk, 
where  k=l,2,...,M,  with  k=M  representing  the  finest  grid,  and  grids  becoming  coarser  as  k 
decreases.  Conventionally,  the  mesh  size  between  two  different  levels  differs  by  a factor  of 
two,  i.e.,  hk+x  = hj 2,  where  hk  represents  the  uniform  grid  spacing  on  Gk-  This  is  not  a 
requirement;  however,  it  does  simplify  the  transfer  of  information  between  adjacent  grid 
levels. 

The  multigrid  algorithm  developed  in  Shyy  and  Sun  (1993)  and  utilized  here,  uses 
the  Full  Multigrid  (FMG)  - Full  Approximation  Storage  (FAS)  concepts  developed  in  Brandt 
(1977)  to  solve  the  governing  system  of  nonlinear  equations.  FAS  refers  to  the  fact  that  at 
each  grid  level  Gk,  the  full  solution  to  the  system  of  equations  is  maintained,  as  opposed  to 
only  the  corrections  to  the  fine  grid  values  on  Gk+i.  This  feature  is  essential  for  solving 
nonlinear  equations  as  it  allows  incremental  changes  in  the  solution  on  coarser  grids  to  be 
reflected  in  both  the  coefficient  matrix  and  source  term  vector  of  the  corresponding  residual 
equation  at  that  level,  which  substantially  increases  the  convergence  rate  of  the  overall 
method.  For  linear  equations,  FAS  is  not  required  since  the  coefficient  matrix  and  source 
term  vector  do  not  depend  on  the  solution,  and  thus,  the  so-called  Correction  Storage  (CS) 
scheme  is  commonly  employed.  To  review  the  CS  strategy,  consider  the  system  of  linear 
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equations  resulting  from  the  discrete  approximation  of  a linear  partial  differential  equation, 
written  as 

AiPk  = fk  (4.1) 

where  Ak  represents  the  coefficient  matrix,  uk  the  exact  solution,  and  fk  the  source  term 
vector  for  multigrid  level  k.  Letting  v*  denote  an  approximate  solution  at  level  k,  such  as 
might  be  obtained  after  a few  relaxation  sweeps  from  an  initial  guess,  the  solution  residual 
rk  is  defined  as 

rk  = Akvk  ~ fk  (4.2) 

Denoting  the  error  vector  as  ek  = uk  - vk,  and  combining  Equation  (4.1)  with  Equation 
(4.2),  an  equation  relating  the  error  and  the  residual  can  be  written  as 

Atfk  = rk  (4-3) 

With  these  expressions,  the  CS  strategy  can  be  summarized  as  follows.  Within  the  first  few 
relaxation  sweeps  at  level  k the  error  vector  ek,  and  thus,  the  residual  vector  rk  are  rapidly 
smoothed;  however,  during  subsequent  sweeps,  little  improvement  is  obtained  since  the  high 
frequency  components  of  the  error  have  already  been  eliminated.  Since  it  is  the  error  and 
residual  which  have  been  smoothed,  the  residual  equation  (not  the  original  equation, 
Equation  (4.1))  is  then  solved  on  the  coarser  level  k— 1.  The  resulting  linear  system  of 
equations  at  level  k-1  is  written  as 

1 = A~'rt  (4.4) 

where  Ak_ ^ represents  the  equivalent  to  matrix  Ak  for  level  k—  1,  and  Ik~  * represents  an 
interpolation  operator  which  converts  a level  k vector  to  a level  k-1  vector.  After  a few 
relaxation  sweeps  on  level  k-1  for  Equation  (4.4),  the  correction  to  the  error  can  be 
transferred  back  to  level  k to  provide  an  improved  estimate  for  the  approximate  solution 


as 
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v 


k 


new 


(4.5) 


The  strategy  illustrated  above  is  effective  for  linear  equations;  however,  it  is 
inadequate  for  nonlinear  equations  since  the  corrections  to  the  solution  on  the  coarser  grid 
levels  cannot  influence  the  coefficient  matrix  and  source  vector,  both  of  which  can  depend 
upon  the  solution  itself.  As  stated  earlier,  the  FAS  strategy  was  devised  to  alleviate  this 
problem.  For  the  FAS  strategy,  the  residual  vector  is  again  defined  by  Equation  (4.2); 
however,  the  coefficient  matrix  Ak  and  the  source  term  vector  fk  are  based  upon  the  solution 
at  the  end  of  the  previous  outer  iteration.  An  outer  iteration  is  here  defined  as  a number  of 
relaxation  sweeps  through  the  linearized  equations  where  the  coefficient  matrix  Ak  and 
source  term  vector  fk  remain  fixed.  Now,  the  discrete  form  of  the  governing  equation  can 
be  written  on  grid  level  k-1  as 


where  the  solution  vector  vJt_1  is  interpreted  as  the  incoming  approximate  solution  from 
level  k plus  the  correction  obtained  by  relaxing  at  level  k-1.  This  is  written  as 


Subtracting  the  restricted  form  of  Equation  (4.2)  from  Equation  (4.6)  results  in  an  equivalent 
nonlinear  form  of  the  residual  equation,  Equation  (4.3)  as 


In  the  current  FAS  strategy,  the  term  Ikk~\Ak_1vk) is  interpreted  as  Ak_  xIk~  xvk where  Ak_  x 
is  computed  once  upon  entering  level  k-1  and  then  held  constant.  Similarly,  Ik~  lfk is  taken 
as  fk_v  which  is  again  computed  only  once  upon  entering  level  k-1.  With  these 
interpretations,  the  final  form  of  the  residual  equation  at  level  k-1  is  written  as 


Bk-\vk-\  ~ Sk- 1 


(4.6) 


vk- 1 = A lyk  + <W-i 


(4.7) 


Bk-i*k-i  - A l(Ak- !»*)  = gk~ i - A~lfk  + A~\  (4.8) 


Bk-Xvk_x  = gk_x  + Ak_x{Ik~ \k)  - fk_x  + (7*-V*) 


(4.9) 
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Since  the  starting  values  for  vk_  j are  l\  1vk , the  source  terms  gk_  j and  fk_  j are  identical 
during  the  first  iteration  at  level  k-1 ; however,  they  differ  during  subsequent  iterations  since 
gk_  j is  continuously  updated  with  the  new  solution  vk_  j . When  the  new  solution  vk_  j has 
been  obtained,  the  solution  at  level  k is  updated  as 

vnew  = vold  + /jfc_i(vn^i  _ (4.10) 

Again,  we  note  that  since  it  is  the  error  components  which  have  been  smoothed  on  level  k, 
it  is  the  changes  in  the  solution  from  relaxation  on  level  k-1  that  are  interpolated  back  to  level 
k,  as  expressed  by  the  second  term  in  Equation  (4.10). 

In  the  current  multigrid  strategy.  Equation  (4.9)  is  used  as  the  residual  equation  on 
the  coarse  grids  for  the  u-momentum  and  v-momentum  equations;  however,  for  the  pressure 
correction  equation,  a different  coarse  grid  residual  equation  is  adopted,  which  is  detailed 
as  follows.  Let  the  pressure  correction  equation  at  multigrid  level  k be  expressed  by  the  linear 
system  of  equations 

AkP\  = bk  (4.11) 

where  Ak  represents  the  coefficient  matrix,  which  is  a function  of  the  control  volume 
dimensions  and  the  velocity  field,  p' k represents  the  pressure  correction  vector,  and  bk 
represents  the  mass  imbalance  vector.  The  corresponding  equation  for  the  coarse  grids  is 
constructed  as  follows.  At  level  k-1,  a pressure  correction  equation  is  still  solved;  however, 
a modified  source  term  is  employed  which  is  composed  of  two  parts,  the  first  being  the 
restricted  mass  imbalance  vector  from  level  k and  the  second  the  change  in  mass  imbalance 
due  to  the  relaxation  of  the  system  of  equations  at  level  k-1.  This  equation  is  written  as 

At-xP\- 1 = lt~'h  + (*?-,  - *1-,)  (4-12) 

The  term  6 ,!_ , represents  the  mass  imbalance  vector  as  computed  from  the  solution  of  the 
momentum  equations  on  level  k-1  during  the  first  outer  iteration  through  the  system  of 


62 


equations.  The  term  b"_  j represents  the  same  quantity,  but  computed  from  the  velocity  field 
obtained  during  the  nth  outer  iteration  at  level  k-1.  Thus,  during  the  first  outer  iteration,  the 
second  term  in  Equation  (4. 1 2)  is  zero,  but  during  subsequent  outer  iterations,  it  contributes 
a change  in  mass  imbalance  to  the  overall  source  term.  Once  the  pressure  field  at  level  k-1 
has  been  determined  from  the  pressure  corrections,  the  pressure  at  level  k is  updated  from 
the  coarser  grid  in  the  same  manner  as  the  velocity  field,  using  an  expression  identical  to  that 
given  in  Equation  (4.10). 

The  FAS  strategy  espoused  here  is  not  a direct  extension  of  that  described  in  Brandt 
(1977);  however,  it  embodies  all  of  the  basic  ideas  associated  with  the  development  of  an 
effective  multigrid  strategy  for  nonlinear  equations.  The  combination  of  Equation  (4.9)  for 
the  coarse  grid  momentum  residual  along  with  Equation  (4.12)  for  the  coarse  grid  pressure 
correction  has  proven  to  be  a robust  technique  for  a variety  of  convection-dominate  flows 
(Shyy  & Sun  1993). 

With  the  discussion  of  FAS  complete,  we  move  now  to  the  Full  Multigrid  (FMG) 
concept.  FMG  refers  to  the  particular  grid  level  cycling  sequence  used  in  progressing  toward 
the  final  converged  solution  on  the  finest  grid.  The  FMG  procedure  implemented  here  is 
based  on  a V cycle,  whereby  one  moves  from  the  finest  grid  to  the  coarsest  grid  and  then  back 
to  the  finest.  A four- level  V cycle  is  shown  below  in  Figure  4.1.  The  FMG  procedure  is 
initiated  on  the  coarsest  grid  level  Gj.  Iterations  are  performed  at  this  level  until  a desired 
level  of  convergence  is  achieved.  This  solution  is  then  used  as  an  initial  guess  for  a two-level 
V cycle  earned  between  Gi  and  G2,  which  proceeds  until  the  desired  convergence  on  G2  has 
been  obtained.  This  solution  is  then  in  turn  used  for  a three-level  V cycle  between  Gi,  G2, 
and  G3.  This  process,  which  is  illustrated  in  Figure  4.2,  is  carried  out  until  convergence  has 
been  reached  on  the  finest  grid  Gm-  In  the  figure,  vk  indicates  the  number  of  relaxation 
sweeps  performed  on  grid  Gk  before  moving  on  to  the  next  level.  In  this  work,  a fixed  number 
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of  relaxation  sweeps  at  each  level  is  performed  for  simplicity;  however,  more  elaborate 
schemes  can  be  implemented  whereby  the  convergence  rate  at  the  current  level  is  monitored 
and  iterations  are  performed  until  converge  becomes  slow,  at  which  time  the  computation 
is  transferred  to  the  next  level. 

In  terms  of  information  transfer  between  the  various  grid  levels  during  the  FMG 
cycling  process,  a combination  of  direct  injection  and  linear  interpolation  is  used  in  this 
work.  With  direct  injection,  if  the  physical  locations  of  any  of  the  dependent  variables 
coincide  for  two  adjacent  grid  levels,  they  are  transferred  directly  between  the  levels  without 
the  use  of  any  interpolation.  If  physical  locations  between  two  adjacent  levels  do  not  coincide 
(which  is  more  often  the  case  due  to  the  nature  of  the  staggered  grid  with  respect  to  grid 
coarsening),  then  linear  interpolation  is  used.  The  use  of  direct  injection  along  with  a 
relatively  low-order  interpolation  for  transferring  data  between  adjacent  grid  levels  is  not 
a critical  issue  in  terms  of  the  basic  effectiveness  of  the  multigrid  technique  because 
effectively  smoothed  information  is  being  transferred,  since  the  high-frequency  components 
have  already  been  eliminated  in  the  relaxation  step.  Regarding  the  transfer  of  the  mass 
imbalance  term  in  the  pressure  correction  equation  from  finer  to  coarser  grids,  a summation 
operator  is  employed,  whereby  the  mass  imbalance  from  multiple  fine  grid  cells  are 


64 


combined  and  used  to  represent  the  incoming  mass  imbalance  on  the  corresponding  coarse 
grid  cell.  For  a complete  discussion  of  interpolation  schemes  for  the  staggered  grid, 
including  the  so-called  full-weighting  scheme,  as  well  as  a complete  development  of  the 
basic  multigrid  methodology  for  the  pressure-based  method  with  the  staggered  grid,  see 
Shyy  and  Sun  (1993)  and  Shyy  et  al.  (1993a). 

4.2  Implementation  Issues  Associated  with  Composite  Multigrid 

In  this  section,  we  examine  two  important  implementation  issues  associated  with  the 
composite  multigrid  solution  technique.  The  first  issue  involves  the  interaction  between  the 
composite  grid  and  multigrid  components.  From  the  discussion,  a methodology  is  chosen 
which  is  considered  most  suitable  for  solving  highly  nonlinear  equations,  such  as  the 
Navier-Stokes  equations,  on  multilevel  composite  grids.  The  second  issue  concerns  the 
coarsening  strategy  which  is  adopted  for  multilevel  overlapping  grids.  This  issue  in 
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considered  from  the  standpoint  of  both  the  numerical  efficiency  achieved  in  obtaining 
solution  to  complex  flows  requiring  composite  grids,  as  well  as  the  flexibility  in  discretizing 
the  domains  for  such  problems. 

4.2.1  Interaction  Between  Composite  Grid  and  Multigrid 

For  multilevel  composite  grids,  one  area  of  primary  interest  involves  the  procedural 
relationship  between  the  composite  grid  component  and  the  multigrid  component.  The 
question  as  to  which  component  should  reside  as  the  outermost  shell  in  solution  algorithms 
for  highly  nonlinear  equations  such  as  the  Navier-Stokes  equations  is  still  largely  unresolved 
and  in  need  of  further  study.  The  two  main  composite  multigrid  approaches  which  have  been 
studied  in  the  literature  can  be  summarized  as  follows.  In  the  first  approach  (hereafter  called 
method  I),  the  composite  grid  component  acts  as  the  outermost  shell.  For  any  block  in  the 
composite  grid,  the  appropriate  boundary  conditions  are  obtained  at  the  finest  grid  level  and 
a multigrid  cycle  is  then  used  to  solve  the  governing  equations  within  the  block.  This 
procedure  is  carried  out  for  each  block  in  the  composite  grid  in  some  predetermined  cyclic 
order.  Thus,  for  each  multigrid  cycle  for  the  composite  grid,  the  boundary  conditions  for 
each  block  are  obtained  only  at  the  finest  grid  level.  In  the  second  approach  (method  II),  the 
multigrid  component  acts  as  the  outermost  shell.  With  this  approach,  for  each  level  in  the 
multigrid  cycle,  the  appropriate  composite  grid  problem  is  solved,  with  internal  boundary 
conditions  for  each  block  being  determined  from  the  neighboring  blocks  within  the  same 
level.  Therefore,  within  each  multigrid  cycle  for  the  composite  grid,  internal  boundary 
information  is  exchanged  multiple  times,  the  exact  number  of  which  depends  on  the  type  of 
cycle  being  executed. 

Model  problem  analysis  was  performed  by  Hinatsu  and  Ferziger  (1991)  to  compare 
the  convergence  rates  obtained  with  the  above  two  methods.  As  test  problems,  they  used  a 


66 


1-D  linear  ODE  and  a 2-D  Poisson  equation.  For  the  1-D  model  problem,  two  subgrids  were 
used  and  both  Dirichlet  and  Neumann  internal  boundary  conditions  (IBCs)  were  employed 
to  transfer  information  between  the  subgrids.  The  influence  of  the  domain  boundary 
conditions  (DBCs)  was  also  considered.  In  some  instances  they  found  that  method  I 
outperformed  method  II  in  convergence  rate;  however,  this  was  highly  dependent  on  both 
the  IBCs  as  well  as  the  DBCs.  In  addition,  the  relative  performance  of  the  two  methods  was 
found  to  be  largely  dependent  on  the  IBC  relaxation  parameter,  which  was  used  to  accelerate 
the  variation  of  the  IBCs.  Overall,  method  II  was  seen  to  be  less  sensitive  to  the  IBC 
relaxation  parameter  and  for  no  relaxation,  slightly  outperformed  method  I.  The  2-D  Poisson 
equation  was  solved  on  a domain  composed  of  two  diagonally  shifted  squares.  Neumann 
conditions  were  imposed  for  both  the  IBCs  and  DBCs.  Overall,  the  fastest  convergence  was 
obtained  using  method  I with  slight  over-relaxation  of  the  IBCs. 

While  method  I was  deemed  more  favorable  than  method  II  in  the  above  study,  others 
have  chosen  to  employ  method  II  for  composite  grids.  Henshaw  and  Chesshire  (1987)  used 
a correction-storage  (CS)  strategy  with  method  II  for  solving  linear,  variable-coefficient 
PDEs.  Tu  and  Fuchs  (1992)  used  method  II  within  a FAS  scheme  based  upon  a V cycle,  for 
solving  the  incompressible  Navier-Stokes  equations  on  staggered  grid  arrangements.  In  their 
work,  multigrid  was  used  for  both  the  momentum  and  pressure  equations,  unlike  other  works 
where  only  the  pressure  equation  is  solved  using  multigrid  (Pemg  & Street  1991).  Their 
results  show  a substantial  speedup  with  the  FAS  scheme,  which  is  largely  maintained  when 
applied  to  composite  grids.  While  the  model  problem  analysis  performed  by  Hinatsu  and 
Ferziger  indicates  a slight  preference  for  method  I,  they  were  employing  a CS  strategy  for 
solving  individual  linear  equations.  Since  we  are  solving  a system  of  nonlinear  equations 
with  a FAS  strategy,  it  is  doubtful  that  their  conclusions  can  be  extended  here. 
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Bearing  the  above  considerations  in  mind,  we  have  chosen  to  employ  method  II  in 
this  work  exclusively.  For  nonlinear  problems,  method  II  allows  changes  in  the  solution 
variables  to  propagate  freely  to  neighboring  blocks  within  the  multigrid  cycle,  thus 
providing  improved  estimates  of  the  boundary  conditions  at  each  grid  level.  This  feature  is 
consistent  with  the  original  FAS  concept  and  allows  the  boundary  values  at  each  grid  level 
to  be  treated  in  the  same  manner  as  the  interior  points.  In  our  opinion,  this  should  enhance 
the  convergence  rate  when  solving  nonlinear  equations  on  composite  grids;  however,  further 
study  is  required  to  determine  if  this  is  indeed  the  case,  and  if  so,  to  what  extent  convergence 
is  enhanced.  Since  the  original  composite  grid  methodology  was  developed  for  general 
systems,  little  difficulty  is  encountered  in  extending  the  original  framework  to  handle 
multilevel  composite  grids.  In  this  regard,  the  original  composite  grid  procedure  operates 
as  a module  within  the  outer  framework  of  the  FMG-FAS  solution  process. 

4.2.2  Coarsening  Strategy  for  Multilevel  Composite  Grids 

Another  issue  which  arises  when  attempting  to  apply  multigrid  to  composite 
overlapping  grids  concerns  the  minimum  overlapping  requirement.  The  minimum 
overlapping  is  determined  by  the  requirement  that  the  sides  of  control  volumes  located  along 
internal  boundaries  must  fall  within  interior  known,  calculated  values  of  the  neighboring 
grid  blocks.  Consider  the  two-block  case  shown  in  Figure  4.3,  where  the  left  and  right  blocks 
have  been  horizontally  separated  for  visualization  purposes,  but  are  assumed  to  overlap 
along  some  vertical  boundary.  Now,  assume  that  the  boundary  conditions  for  block  2 are  to 
be  determined.  In  this  case  the  required  overlap  is  determined  only  by  the  left  boundary  for 
block  2 and  the  right  boundary  for  block  1 , since  the  surrounding  boundaries  are  Dirichlet 
conditions.  For  a vertical  internal  boundary  between  the  adjacent  blocks,  the  two  types  of 
boundary  control  volumes  of  interest  are  the  continuity  and  v-velocity  component  control 
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volumes,  as  detailed  in  Chapter  3.  These  two  types  of  control  volumes  along  the  left  internal 
boundary  of  block  2 are  shown  in  the  figure.  For  the  continuity  control  volumes  along  the 
left  internal  boundary  of  block  2,  the  inlet  u velocity  component  profile  along  the  line 
denoted  by  A2  is  the  required  boundary  condition.  This  line  must  fall  within  interior, 
calculated  u velocity  component  values  of  block  1 and  with  the  staggered  grid  convention 
means  that  line  A2  must  fall  on  or  to  the  left  of  line  Aj  in  block  1.  Boundary  conditions  for 
the  v control  volumes  require  the  calculation  of  the  v-momentum  flux  from  block  1 which 
requires  both  u-  and  v-  velocity  component  values.  Now,  the  rightmost  interior,  calculated 
v-component  values  in  block  1 are  located  along  line  Bv  Since  u-component  values  are  also 
required,  the  minimum  overlapping  again  requires  that  line  A2  of  block  2 fall  on  or  to  the 
left  of  line  A j in  block  1 , identical  to  the  case  of  the  control  volume  used  to  compute  the 
continuity  equation.  Thus,  based  on  these  requirements,  line  Al  in  block  1 defines  the 
minimum  overlapping  requirement  for  block  2.  Similarly,  line  B2  defines  the  minimum 


Block  1 Block  2 


Figure  4.3.  Two-block  composite  grid  used  to  demonstrate  minimum 
overlapping  requirement. 
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overlapping  requirement  for  block  1 . For  a valid  overlap,  both  overlapping  constraints  must 
be  satisfied.  Similar  constraints  also  exist  for  horizontal  overlapping  zones. 

With  the  multigrid  procedure  adopted  in  Shyy  and  Sun  (1993),  successively  coarser 
grids  are  created  by  removing  every  other  interior  grid  line  in  both  coordinate  directions 
from  the  grid  at  the  previous  (finer)  level.  Letting  M denote  the  total  number  of  multigrid 
levels  and  nk  the  number  of  grid  lines  in  any  coordinate  direction  for  grid  level  k,  then 
nk  = (nk+ 1 + l)/2for  1 < k < M — 1.  With  this  notation,  level  Mis  again  the  finest  grid 
and  coarser  grids  have  lower  k values.  For  composite  overlapping  grids,  this  type  of 
coarsening  strategy  quickly  leads  to  a violation  of  the  overlapping  requirements  just 
discussed  unless  appropriate  measures  are  taken.  Consider  a two-block  continuous  grid  line 
case  as  shown  below  in  Figure  4.4(a).  The  two  blocks  combine  to  form  a vertical  single  cell 
overlap  region  which  satisfies  the  overlapping  requirements  for  the  finest  level.  However, 


(a)  (b) 


Figure  4.4.  Two-block  composite  overlapping  grid  used  to  demonstrate 
violation  of  overlap  constraint  with  grid  coarsening,  (a)  Valid 
overlap  at  finest  level;  (b)  Invalid  overlap  after  coarsening. 
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at  the  next  multigrid  level  shown  in  Figure  4.4(b),  the  overlapping  requirements  are  not  met, 
and  thus,  the  overlap  is  not  valid.  In  general,  for  a composite  grid  with  uniform  grid  spacing 
and  continuous  grid  lines  across  internal  boundaries,  in  order  to  ensure  a valid  overlap  on 
the  coarsest  multigrid  level,  an  overlap  thickness  of  2M~ 1 cells  must  be  used  on  the  finest 
grid  level.  In  the  FMG  procedure,  for  M=6,  an  overlap  thickness  of  32  cells  is  required. 
Clearly,  this  strategy  is  very  wasteful  of  computational  time  and  physical  memory,  as  many 
redundant  grid  nodes  are  required  just  to  maintain  the  appropriate  overlap. 

One  possible  remedy  to  this  situation  hinges  on  the  fact  that  with  the  above  strategy, 
the  location  of  the  physical  boundaries  of  each  grid  block  were  maintained  during  the 
coarsening  process.  For  some  simple  composite  grid  configurations,  such  as  a two-block 
cavity  flow  problem,  it  is  possible  to  adjust  the  physical  location  of  certain  internal  block 
boundaries  at  each  multigrid  level  to  ensure  that  the  appropriate  overlap  is  maintained.  Thus 
for  each  multigrid  level,  the  physical  locations  of  the  internal  boundaries  are  different. 
Algorithmically,  this  type  of  strategy  is  not  difficult  to  implement  for  simple  configurations 
such  as  the  one  above;  however,  for  configurations  in  which  multiple  ( >3)  blocks  overlap 
at  the  same  location,  or  in  which  some  physical  boundaries  are  composed  of  both  boundary 
condition  segments  and  internal  boundary  segments,  this  strategy  becomes  very  difficult  to 
implement  and  in  some  cases  is  impossible  due  to  the  fact  that  physical  boundary  constraints 
may  be  violated  by  moving  sides  with  internal  boundaries. 

A general  solution  to  this  problem  lies  in  modifying  the  coarsening  method  used 
during  the  restriction  process.  With  the  previously  mentioned  coarsening  method,  the  grid 
at  level  k-1  is  created  from  level  k by  removing  every  other  interior  grid  line  in  both 
coordinate  directions.  As  a result  of  this  procedure,  the  rows  of  cells  directly  adjacent  to  the 
boundaries  of  the  grid  thicken  as  k -»  1 until  eventually  the  off-boundary  side  of  the  row 
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leaves  the  overlap  region  established  by  the  finest  grid  level,  thus  creating  a non-valid 
overlap.  In  order  to  prevent  this,  a row  or  rows  of  fine  grid  cells  can  be  retained  next  to  the 
boundaries  at  each  multigrid  level.  For  composite  grid  schemes  utilizing  standard  5-point 
stencils  (such  as  the  hybrid  scheme  for  the  convection  terms  and  the  central  difference 
scheme  for  diffusion  terms),  it  is  sufficient  to  maintain  a single  layer  of  fine  grid  cells  at  the 
boundaries  for  all  multigrid  levels.  For  9-point  stencils  commonly  arising  from  the  use  of 
higher-order  convection  schemes,  at  least  two  layers  of  fine  grid  cells  must  be  maintained. 
In  this  work,  we  will  retain  two  rows  of  fine  grid  cells  at  each  multigrid  level.  With  this 
coarsening  strategy,  when  employing  5-point  stencil  schemes,  now  for  a composite  grid  with 
uniform  grid  spacing  and  continuous  grid  lines  across  internal  boundaries,  a single-cell 
overlap  thickness  on  the  finest  grid  level  is  sufficient  to  ensure  a valid  overlap  at  all  multigrid 
levels,  while  for  9-point  stencils  a two-cell  overlap  is  sufficient.  Letting  n^  again  denote 
the  number  of  grid  lines  in  any  coordinate  direction  for  multigrid  level  k,  then 
nk  = ( nk+l  + 5)/2  for  1 < k < M — 1.  As  an  example  of  this  grid  coarsening 
strategy,  consider  the  four- level  multigrid  for  a square  lid-driven  cavity  shown  in  Figure  4.5. 
The  grid  shown  in  Figure  4.5(a)  is  the  finest  grid,  and  Figures  4.5(b)-(d)  represent  successive 
coarsenings  from  previous  levels.  The  finest  grid  contains  21x21  uniformly-spaced  nodes, 
while  successive  grids  contain  13x13,  9x9,  and  7x7  nodes  respectively.  It  should  be  noted 
that  this  type  of  uneven  grid  coarsening  procedure  has  been  developed  previously  in  a single 
block  multigrid  procedure  (Shyy  etal.  1993a).  The  primary  motivation  there  was  to  maintain 
adequate  grid  resolution  in  the  region  close  to  the  solid  boundary.  For  high  Reynolds  number 
flows,  especially  those  in  the  turbulent  regime,  it  has  been  found  that  this  type  of  grid 
coarsening  procedure  can  help  maintain  the  necessary  length  scale  resolution  requirement, 
resulting  in  improved  convergence  rates  with  the  multigrid  solver.  Here,  this  aspect  will  still 


be  useful. 
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4.3  Evaluation  of  Composite  Multigrid  Technique 

In  this  section,  three  test  cases  are  used  to  examine  the  effectiveness  of  the  composite 
multigrid  technique  which  has  been  developed.  First  we  consider  a simple  geometry  in  order 
to  investigate  the  interaction  between  the  composite  gridding  procedure  and  the  multigrid 
technique.  The  well-known  lid-driven  cavity  flow  is  used  as  the  test  case.  For  the  next  two 
cases,  we  compute  two  complex  flows,  one  the  baffle  flow  previously  presented  in  Chapter 


(c)  (d) 


Figure  4.5.  Four-level  multigrid  for  square  lid-driven  cavity  showing  new 
coarsening  strategy  for  composite  multigrid,  (a)  Level  four; 
(b)  Level  three;  (c)  Level  two;  (d)  Level  one. 
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3 and  the  other  the  flow  about  a CFD  configuration.  These  cases  are  used  to  assess  the  ability 
of  the  current  procedure  to  efficiently  compute  fluid  flows  with  complex  geometry. 

4.3.1  Influence  of  Composite  Gridding  on  Multigrid  Performance 

For  composite  grid  methods,  two  factors,  namely,  the  convergence  rate  within  the 
individual  subgrids  and  the  presence  of  the  internal  boundaries,  interact  to  determine  the  net 
effectiveness  of  the  method.  For  a fixed  physical  domain  and  grid  resolution,  as  the  domain 
is  divided  into  a larger  number  of  smaller  subgrids,  the  convergence  rate  within  each  of  the 
subgrids  increases.  At  the  same  time,  however,  the  increased  presence  of  the  internal 
boundaries  reduces  the  coupling  between  the  unknowns  in  the  different  subgrids  and 
eliminates  the  implicit  connection  to  boundary  conditions  for  some  subgrids.  Thus,  while 
the  convergence  rate  of  the  linearized  equations  within  individual  subgrids  may  be 
improved,  that  of  the  overall  system  may  not.  The  net  effect  of  composite  gridding  is 
certainly  problem  dependent;  however,  for  the  current  procedure,  it  has  been  shown  (Shyy 
et  al.  1993b)  that  for  a square,  lid-driven  cavity  flow,  the  convergence  rate  is  virtually 
unaffected  by  the  way  in  which  the  domain  is  partitioned,  or  the  number  of  subgrids  used. 
At  the  same  time,  however,  for  a high  aspect  ratio  channel  flow,  division  of  the  domain  into 
several  subgrids  was  found  to  significantly  increase  the  rate  of  convergence  of  the  overall 
calculation  (by  almost  an  order  of  magnitude  for  large  aspect  ratio). 

With  the  incorporation  of  the  composite  grid  method  into  the  outer  framework  of  a 
multigrid  strategy,  the  interaction  of  the  composite  grid  and  multigrid  components  becomes 
a relevant  issue  of  concern.  To  examine  this  interaction  and  the  overall  effect  of  composite 
gridding  on  multigrid  efficiency,  a single-block,  lid-driven  cavity  flow  is  computed  and 
compared  with  corresponding  multi-block  configurations  for  the  same  flow.  The  single-  and 
multi-block  grids  used  are  shown  in  Figure  4.6.  Grid  A consists  of  a single  block  of  69x69 
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nodes,  while  grids  B and  C consist  of  two  69x37  blocks  and  four  69x21  blocks  respectively. 
These  grid  resolutions  were  chosen,  because  with  a 69x69  initial  multigrid  level,  five 
coarsenings  can  be  achieved  with  the  new  coarsening  procedure,  resulting  in  successive 
grids  of  37x37,  21x21,  13x13,  9x9,  and  7x7  nodes.  For  cases  B and  C,  the  resulting 
configurations  have  overlaps  of  five  cells  and  four  cells  respectively  on  the  finest  grid  level; 
however,  as  stated  previously,  for  this  case  the  convergence  rate  is  essentially  unaffected  by 
the  degree  of  overlap,  and  thus  any  overlap  satisfying  the  minimum  overlapping 
requirements  is  sufficient  for  this  case. 

Since  the  multi-block  grids  contain  more  control  volumes  than  the  single  grid, 
equivalent  convergence  tolerances  must  be  used.  For  the  following  computations,  a 
normalized  residual  of  1.0  X 10  ~4  is  used  for  the  single  grid  case.  The  corresponding 
normalized  residuals  for  cases  B and  C are  1.07  x 10 ~4  and  1.21  x 10  ~4  respectively. 
These  normalized  global  residuals  are  used  to  give  and  indication  of  when  the  various  cases 
have  reached  the  same  level  of  convergence.  Another  note  needs  to  be  made  regarding  the 


A 


69x69 


B 

C 

69x37 

69x21 

69x21 

69x37 

69x21 

69x21 

Figure  4.6.  Grids  used  for  testing  effect  of  composite  gridding  on 
multigrid  performance. 
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relative  amount  of  work  performed  between  successive  multigrid  levels.  In  the  following 
results,  the  amount  of  work  required  to  achieve  convergence  is  stated  in  work  units.  A work 
unit  is  here  defined  as  one  sweep  through  all  of  the  blocks  on  the  finest  grid.  Since  we  are 
not  using  the  standard  coarsening  procedure,  whereby  the  number  of  grid  points  at  any 
multigrid  level  k is  half  that  at  level  k+1,  the  work  per  iteration  on  grid  level  k is  not  one 
fourth  of  that  on  level  k+1  per  iteration.  For  the  new  coarsening  procedure,  the  work  required 
per  iteration  on  level  k is  greater  than  one  fourth  of  that  on  level  k+1  and  is  dependent  on 
the  configuration.  For  the  three  cases  to  be  considered  here,  the  work  per  iteration  on  each 
of  the  possible  multigrid  levels  with  respect  to  the  finest  grid  level  is  shown  in  Table  4.1. 


Table  4.1.  Relative  work  for  successive  multigrid  levels  compared  to  finest  grid  level 
for  three  cavity  flow  cases.  Note:  Only  five  multigrid  levels  are  possible 
for  case  B and  four  for  case  C. 


Level 

Number 

Relative  Work  Compared  to  Finest  Level 

Case  A 

Case  B 

Case  C 

6 

1.0000 

5 

0.2875 

1.0000 

4 

0.0926 

0.3043 

1.0000 

3 

0.0355 

01069 

0.3320 

2 

0.0170 

0.0458 

0.1304 

1 

0.0103 

0.0247 

0.0628 

In  the  FMG  cycle,  the  number  of  relaxation  sweeps  at  each  level  are  taken  as  follows. 
For  case  A,  (vfe)fc= 16  = 12, 8, 8, 4, 4, 2 . Thus,  for  a two-level  cycle  for  case  A,  levels  five 
and  six  are  active,  with  v5  = 4 and  v6  = 2.  The  corresponding  values  for  cases  B and  C 
are  respectively,  (v*)*=15  = 8, 8, 4, 4, 2 and  (vk)k=14  = 8,4,4, 2.  During  the  course  of 
the  FMG  cycle  the  current  level  is  taken  to  complete  convergence,  i.e.  the  equivalent 
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convergence  tolerance  of  1 .0  x 10  ~4,  before  the  next  multigrid  level  is  initiated.  While  this 
requirement  is  wasteful  of  iterations  on  coarser  grid  levels  during  the  initial  process  of 
obtaining  an  approximate  solution  for  the  finest  grid  V cycle,  it  does  not  affect  the  results 

to  be  presented. 

Results  for  cases  A,  B,  and  C are  presented  in  Tables  4.2.  For  each  case,  the  number 
of  work  units  required  to  achieve  convergence  and  the  corresponding  speedup  with  respect 
to  the  single  level  calculation,  in  terms  of  a work  unit  ratio,  are  tabulated  for  each  possible 
number  of  multigrid  levels.  Reynolds  numbers  of  both  100  and  1000  were  computed  in  order 
to  evaluate  the  flow  dependency  of  the  results  as  well  as  to  provide  a database  for 
computations  in  the  following  section  for  complex  flows.  In  all  of  the  calculations,  the 
hybrid  scheme  (Spalding  1972)  was  used  for  the  convection  terms  while  standard  central 
differences  were  employed  for  the  diffusion  terms.  The  hybrid  scheme  was  chosen  in  order 
to  have  the  same  scheme  at  all  multigrid  levels,  since  with  the  central  difference  scheme, 
divergence  of  the  solution  was  found  to  occur  at  some  coarser  grid  levels  for  Re=1000.  In 
any  case,  the  hybrid  scheme  defaults  to  the  central  difference  scheme  for  a majority  of  the 
flow  field  for  both  Re=100  and  Re=1000  since  the  local  cell  Peclet  number  exceeds  the  value 
of  2 on  the  finest  grid  only  near  the  upper  corners  of  the  domain  (Thakur  1993). 

From  the  tabulated  values  shown  in  Tables  4.2,  several  observations  can  be  made. 
Firstly,  we  see  that  the  number  of  work  units  required  to  achieve  convergence  for  the 
single-level  computations  for  cases  A,  B,  and  C are  virtually  the  same,  for  both  Re=100  and 
Re=1000,  indicating  the  relative  insensitivity  of  the  computations  to  the  number  of  blocks 
used  to  partition  the  domain.  This  is  consistent  with  previous  findings  of  Shyy  et  al.  (1993b) 
for  the  same  flow,  but  with  different  grid  resolutions  and  domain  partitioning.  Secondly,  for 
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Tables  4.2.  Multigrid  computation  results  for  cavity  flow, 
(a)  Case  A;  (b)  Case  B;  (c)  Case  C. 

(a) 


Re  = 

100 

Re  = 

1000 

MG  Levels 

Work  Units 

Speedup 

Work  Units 

Speedup 

1 

4978 

1.00 

1950 

1.00 

2 

1022 

4.87 

1527 

1.28 

3 

616 

8.08 

1253 

1.56 

4 

524 

9.50 

1113 

1.75 

5 

548 

9.09 

1145 

1.70 

6 

580 

8.58 

1203 

1.62 

(b) 


Re  = 100 

Re  = 1000 

MG  Levels 

Work  Units 

Speedup 

Work  Units 

Speedup 

1 

4976 

1.00 

2099 

1.00 

2 

1076 

4.62 

1611 

1.30 

3 

673 

7.40 

1335 

1.57 

4 

596 

8.35 

1197 

1.75 

5 

618 

8.06 

1164 

1.80 

(c) 


Re  = 100 

Re  = 1000 

MG  Levels 

Work  Units 

Speedup 

Work  Units 

Speedup 

1 

5030 

1.0 

2152 

1.00 

2 

1188 

4.24 

1696 

1.27 

3 

911 

5.52 

1449 

1.49 

4 

824 

6.10 

oscillatory 

this  configuration,  the  solution  for  Re=1000  converges  much  quicker  than  that  for  Re=100. 
This  surprising  result  was  also  observed  by  Shyy  and  Sun  (1993)  and  is  inconsistent  with 
other  configurations  whereby  the  number  of  iterations  increases  as  the  Reynolds  number  is 


increased. 
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Regarding  the  effect  of  composite  gridding  on  multigrid  performance,  we  see  that 
as  the  number  of  subgrids  increases,  multigrid  effectiveness  decreases.  For  two  multigrid 
levels,  the  reduction  in  effectiveness  is  only  slight;  however,  as  the  number  of  multigrid 
levels  increases,  the  effectiveness  is  further  reduced.  This  can  be  explained  as  follows.  For 
the  single  grid  case  A,  the  finest  grid  level  contains  69x69  nodes,  resulting  in  successive 
multigrid  levels  of  37x37,  21x21,  13x13,  9x9,  and  7x7  nodes.  As  the  number  of  multignd 
levels  is  increased  for  this  case,  we  see  a continual  improvement  in  convergence  rate  up  to 
five  multigrid  levels,  at  which  point  the  convergence  rate  decreases.  This  clearly  results  from 
the  relative  resolutions  of  multigrid  levels  three,  two,  and  one,  which  are  13x13,  9x9,  and 
7x7,  respectively.  Since  these  resolutions  are  comparable,  the  errors  introduced  via 
interpolation  during  the  prolongation  and  restriction  process  are  not  efficiently  eliminated 
since  the  wavelengths  of  the  errors  are  comparable  to  the  grid  spacing  onto  which  they  are 
being  interpolated.  In  other  words,  the  additional  speedup  obtained  with  the  introduction  of 
a new  multigrid  level  will  be  little  or  none,  if  the  resolution  of  the  new  level  is  comparable 
to  that  of  the  previous  level.  For  case  B,  which  starts  with  an  initial  subgrid  resolution  of 
69x37,  this  limitation  is  reached  with  fewer  multigrid  levels  than  case  A,  due  to  the 
resolution  of  only  37  nodes  in  the  vertical  direction  on  each  of  the  finest  subgrids.  Case  C, 
which  starts  with  an  initial  subgrid  resolution  of  69x21,  approaches  this  limitation  even 
sooner.  Thus,  for  the  three  grid  configurations  A,  B,  and  C,  the  difference  in  resolution 
between  any  two  multigrid  levels  is  the  least  for  case  C,  followed  by  case  B,  and  then  case 
A.  Therefore,  a reduction  in  multigrid  effectiveness  should  occur  as  the  domain  is  partitioned 
into  two  blocks,  as  in  case  B,  and  then  four  blocks,  as  in  case  C. 

The  situation  described  above,  in  which  newly  introduced  multigrid  levels  have 
roughly  the  same  resolution  as  their  predecessors,  is  a direct  consequence  of  the  new 
coarsening  strategy,  and  introduces  a tradeoff  between  domain  partitioning  and  the 
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computational  efficiency  which  can  be  obtained  with  composite  multigrid.  For  the  cases 
computed  above,  however,  the  grid  resolutions  employed  on  the  finest  levels  are  still 
relatively  coarse,  and  thus  the  ratio  of  successive  multigrid  level  resolutions  quickly 
decreases  below  a factor  of  two  (which  is  always  maintained  for  the  standard  coarsening 
procedure)  as  the  number  of  multigrid  levels  is  increased.  For  much  finer  grid  resolution 
cases,  for  which  the  multigrid  technique  is  most  useful,  little  difference  exists  between  the 
new  coarsening  procedure  and  the  standard  procedure,  and  the  tradeoff  between  domain 
partitioning  and  multigrid  efficiency  is  minimal.  In  terms  of  this  tradeoff,  one  last 
observation  from  the  cases  presented  above  is  noted.  While  the  overall  multigrid 
effectiveness  is  much  less  for  all  cases  with  Re=1000,  the  sensitivity  of  the  effectiveness  to 
domain  partitioning  is  less  for  Re=1000  than  for  Re=100.  For  two  multigrid  levels  the 
degradation  in  multigrid  performance  between  cases  A and  C is  13%  for  Re=100,  and  1% 
for  Re=1000,  while  for  three  multigrid  levels  a degradation  of  32%  occurs  for  Re=100,  and 
only  5%  occurs  for  Re=1000. 

4.3.2  Composite  Multigrid  for  Complex  Flows 

In  order  to  test  the  current  composite  multigrid  procedure  for  complex  flows,  two 
multiply  connected  flows  are  computed  here.  The  first  case  considered  is  the  baffle  flow 
previously  considered  in  Chapter  3.  The  flow  is  computed  for  Reynolds  numbers  of  100  and 
800  based  on  the  inlet  velocity  and  configuration  width.  The  hybrid  scheme  of  Spalding 
(1972)  is  employed  here  for  both  Reynolds  numbers.  The  finest  multigrid  level  corresponds 
identically  to  the  case  computed  in  Chapter  3.  Due  to  the  x-coordinate  resolution  of  41  nodes 
in  block  four,  at  most  three  multigrid  levels  can  be  used  for  the  computation.  Figure  4.7 
shows  the  finest  grid  level  along  with  the  two  successive  multigrid  coarsenings  obtained 
using  the  new  coarsening  strategy.  We  note  here  that  grid  lines  within  some  blocks  of  the 
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coarsest  multigrid  level  have  been  adjusted  so  as  to  guarantee  global  mass  conservation  for 
each  of  the  those  blocks,  as  discussed  in  Section  3.5. 

Two  important  comments  must  be  made  in  terms  of  the  multigrid  implementation  for 
this  case.  Regarding  the  pressure  cycling  procedure  mentioned  in  Chapter  3 for  ensuring  a 
continuous  pressure  field,  in  order  to  obtain  full  convergence  it  was  found  that  this  procedure 
should  only  be  performed  on  the  finest  multigrid  level.  Application  of  pressure  cycling  at 
all  multigrid  levels  produced  an  oscillatory  convergence  path  which  never  achieved  full 
convergence.  No  concrete  explanation  for  this  behavior  can  be  given  at  the  current  time.  The 
other  comment  pertains  to  the  coarse  grid  pressure  correction  equation,  Equation  (4.12). 
With  the  pressure  cycling  strategy,  the  pressure  correction  residuals  (mass  imbalances) 
within  block  six,  at  the  interfaces  with  blocks  four  and  five  were  found  to  be  roughly  an  order 
of  magnitude  higher  than  those  throughout  the  rest  of  the  domain  during  the  course  of  the 
entire  calculation.  This  is  caused  by  the  introduction  of  the  pressure  discontinuities  along 
those  lines  within  block  six  during  the  intermediate  iterations  before  convergence.  With  the 
application  of  Equation  (4.12)  for  the  coarse  grid  pressure  correction  these  residuals 
propagate  throughout  the  entire  multigrid  cycle  and  severely  hinder  the  effectiveness  of  the 
procedure.  In  order  to  alleviate  this  problem,  the  mass  imbalances  between  multigrid  levels 
were  effectively  decoupled  by  removing  the  first  term  from  the  right-hand  side  of  Equation 
(4.12)  for  the  coarse  grid  pressure  corrections.  This  resulted  in  a substantial  improvement 
in  multigrid  effectiveness  (by  a factor  of  two  in  terms  of  speedup  ratio  for  this  case),  and  thus, 
for  all  of  the  following  results  this  new  procedure  has  been  employed.  It  is  interesting  to  note 
that  for  a single-block,  lid-driven  cavity,  virtually  no  difference  in  speedup  ratio  occurs 
between  the  original  method  and  this  new  method. 

Results  from  the  multigrid  computations  for  both  Re= 1 00  and  Re=800  are  presented 
in  Table  4.3,  in  terms  of  the  number  of  fine  grid  iterations  required  to  achieve  convergence 
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on  the  finest  grid  level,  the  total  work  performed  during  the  FMG  cycle  and  the  speedup  in 
terms  of  a work  unit  ratio.  Figure  4.8  shows  a comparison  of  the  convergence  paths  of  the 
single-level  and  multi-level  computations  for  both  Reynolds  numbers.  The  speedup  ratios 
obtained  for  this  configuration  demonstrate  that  even  for  complex  composite  grid 
configurations,  multigrid  effectiveness  can  be  largely  maintained.  In  general,  for  any 
composite  grid  configuration,  the  maximum  attainable  speedup  ratio  will  be  determined 
roughly  by  the  lowest  resolution  block  of  the  composite  grid.  For  this  configuration,  the 
lowest  resolution  grid  blocks  are  blocks  one  and  seven,  each  having  a resolution  of  21x53 
nodes.  Therefore,  one  would  expect  the  speedup  ratio  to  fall  close  to  that  for  a 21x21  single 
block  cavity  flow,  which  is  consistent  with  the  findings  of  Shyy  and  Sun  (1993).  Finally,  the 
multigrid  effectiveness  obtained  for  this  configuration  is  again  dependent  upon  Reynolds 
number  and  follows  the  same  trends  observed  for  the  single-block,  lid-driven  cavity  flow 
previously  presented. 


Table  4.3.  Multigrid  computations  for  baffle  flow. 


Re  = 100 

Re  = 800 

# levels 

# f.  g.  iter. 

# work  u. 

speedup 

# f.  g.  iter. 

# work  u. 

speedup 

1 

1977 

1977 

1.00 

1288 

1288 

1.00 

2 

720 

1099 

1.80 

660 

1010 

1.28 

3 

242 

612 

3.23 

340 

835 

1.54 

Now  we  consider  a flow  domain  composed  of  a cavity  with  the  internal  obstacle 
CFD.  The  physical  configuration  with  boundary  conditions  is  shown  in  Figure  4.9(a).  The 
grid  has  been  constructed  in  such  a manner  as  to  provide  not  only  a sufficient  number  of 
multigrid  levels  for  each  subgrid,  but  subgrids  which  do  not  violate  the  global  mass 
conservation  property  upon  coarsening.  Thirty  one  subgrids  are  used  to  discretize  the 
domain.  Figure  4.9(b)  shows  the  subgrid  partitioning  used.  Overlap  regions  have  not  been 
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shown  for  clarity.  Each  of  the  subgrids  is  composed  of  either  21x21, 21x37, 37x21,  or  37x37 
nodes,  so  that  four  multigrid  levels  are  possible  for  this  configuration.  The  different 
multigrid  levels  are  shown  in  Figure  4. 1 0. 

The  flow  in  the  main  portion  of  the  enclosure  is  computed  for  Re=100  based  on  the 
sliding  lid  velocity  and  configuration  width.  For  the  enclosed  lid  driven  cavity  in  the  letter 
D,  the  lid  velocity  is  the  same  as  the  main  configuration,  resulting  in  Re=6  based  on  the 
enclosed  cavity  width.  Pressure  cycling  is  again  required  at  every  outer  iteration  and  is 
performed  in  a similar  manner  to  that  described  in  Chapter  3.  As  with  the  baffle  flow, 
pressure  cycling  is  carried  out  only  on  the  finest  grid  level.  Streamlines  of  the  flow  computed 
on  the  finest  grid  level  are  shown  in  Figure  4.11.  Results  for  the  various  multigrid 
computations  are  presented  in  Table  4.4  in  terms  of  the  number  of  fine  grid  iterations 
required  to  reach  convergence,  the  corresponding  number  of  work  units  for  the  entire  FMG 
cycle  and  the  effective  speedup  in  terms  of  a total  work  unit  ratio.  Since  the  resolution  for 
this  configuration,  based  on  the  total  number  of  grid  nodes,  is  equivalent  toal71xl71  square 
grid,  convergence  was  taken  to  a normalized  tolerance  of  10  ~3  . From  the  speedup  ratio 
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u = 1 


u=  1 


(b) 


Figure  4.9.  CFD  configuration,  (a)  Physical  configuration  with  boundary 
conditions;  (b)  Subgrid  partitioning  for  multigrid  CFD 
configuration. 


values  shown  in  Table  4.4,  we  see  that  the  multigrid  procedure  is  still  very  effective  in 
reducing  the  total  amount  of  work  required  to  achieve  convergence  on  the  finest  grid  level. 
In  fact,  a larger  speedup  ratio  is  obtained  for  this  configuration  than  for  the  baffle  flow,  due 
primarily  to  the  finer  discretization  on  the  finest  grid  level.  These  results  indicate  that  for 
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(a) 


(b) 


Figure  4.10.  Multigrid  levels  for  CFD  configuration,  (a)  Level 
one;  (b)  Level  two;  (c)  Level  three;  (d)  Level  four. 
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(c) 


(d) 


Figure  4.10  — continued 
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high-resolution  grids  (much  higher  than  computed  here)  we  can  anticipate  an  order  of 
magnitude  speedup  even  for  very  complex  multiply-connected  configurations. 


Table  4.4.  Multigrid  computations  for  CFD  configuration. 


# levels 

# fine  grid  it. 

# work  units 

speedup 

1 

2735 

2735 

1.00 

2 

1020 

1740 

1.57 

3 

468 

1284 

2.13 

4 

228 

782 

3.50 

4.4  Concluding  Remarks 

In  this  chapter,  several  issues  associated  with  the  implementation  of  a composite 
multigrid  procedure  for  the  pressure-based  solution  of  the  Navier-Stokes  equations  have 
been  addressed.  A FAS  strategy  is  adopted  here  and  implemented  into  a FMG  framework 
whereby  the  composite  grid  procedure  acts  as  a module  within  an  outer  multigrid  shell.  A 
new  coarsening  strategy  has  been  presented,  which  is  amenable  for  composite  overlapping 
grids  in  terms  of  both  the  flexibility  in  discretization  as  well  as  solution  efficiency.  Using  a 
simple  lid-driven  cavity  flow  problem,  it  was  shown  that  some  degradation  in  multigrid 
efficiency  can  occur  with  the  new  coarsening  strategy  for  composite  grids  where  partitioning 
of  the  domain  results  in  relatively  coarse  grids  at  the  finest  grid  level.  However,  as  the  overall 
grid  resolution  is  increased,  for  which  the  multigrid  method  is  most  effective,  the  new 
coarsening  strategy  approaches  the  standard  method  most  commonly  employed  and  the 
tradeoff  between  domain  partitioning  and  multigrid  effectiveness  disappears.  Two  complex 
flows  were  computed  with  the  composite  multigrid  procedure  and  used  to  demonstrate  the 
effectiveness  of  the  overall  method.  Results  for  these  flows  indicate  that  the  multigrid 
effectiveness  is  largely  maintained,  and  thus,  the  current  procedure  is  shown  to  provide  a 
flexible  and  efficient  means  of  computing  truly  complex  fluid  flows. 


CHAPTER  5 

COMPOSITE  GRID  METHOD  FOR  SCALAR  TRANSPORT 


5.1  Introduction 

In  the  previous  chapters,  a composite  grid  methodology  has  been  developed  for 
isothermal  flows  governed  by  the  incompressible  form  of  the  Navier-Stokes  equations.  In 
developing  a composite  grid  method  for  flows  of  variable  temperature,  or  flows  involving 
additional  scalar  quantities,  new  issues  arise,  primarily  associated  with  the  treatment  of  the 
internal  boundaries  for  the  auxiliary  equations  which  must  be  solved  (e.g.,  energy  equation). 
Here,  we  consider  an  extension  of  the  existing  methodology  to  the  treatment  of  flows  driven 
by  natural  convection,  for  which  the  auxiliary  equations  are  those  governing  field  properties 
which  directly  affect  the  density.  Such  flows  encompass  many  of  the  issues  which  arise  in 
the  computation  of  more  complicated  flows  and  can  serve  as  a basis  for  discussion  for  several 
key  points  regarding  the  development  of  composite  grid  techniques. 

Basic  analysis  and  development  of  a composite  grid  methodology  for  natural 
convection  flows  will  be  carried  out  for  flows  in  which  properties  such  as  temperature  and 
salinity  contribute  to  the  density  variation.  The  governing  equations  of  interest  here  are  the 
time-dependent  forms  of  the  continuity,  momentum,  energy  and  salinity  equations  obtained 
after  employing  the  Boussinesq  approximation.  These  can  be  written  as 

+^(0ov)  = 0 (5.1a) 
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jjteo") + £<00““)  + £<0o“v)  - - £ + £<"£> + £<“§) 

ft<e0v>  + £(0o“v)  + £(6,ovv)  = - aj  + + ty&p  + B 

f + £w  + £w  = £(«rf)  + £(arf) 

f + £w> + £w  - £<a*f> + £<a*f > 1 

where  ar  is  the  thermal  diffusivity  and  as  is  the  mass  diffusivity  of  the  salt.  The  term  B 
in  the  vertical  momentum  equation,  Equation  (5.  lc),  is  the  source  term  due  to  buoyancy,  and 
is  given  by  the  expression 

B = QoSWAT  ~ To)  ~ Ps(S  ~ S0)]  (5.2) 

The  Boussinesq  approximation  assumes  that  the  fluid  density  g,  varies  only  slightly  from 
a reference  density  state  ^0>  a°d  can  be  expressed  as 

P = Po  - QW  ~ T0)  - faS  - So)]  (5.3) 

In  addition,  the  density  variations  are  only  assumed  to  produce  the  fluid  motion  via  the 
source  term  B,  and  are  not  considered  to  contribute  to  any  of  the  other  terms  in  the  equations. 
These  assumptions  are  valid  when  the  resulting  temperature  and  salinity  variations  do  not 
produce  large  density  variations.  For  the  cases  considered  here,  variations  in  the  density  of 
less  than  1%  are  expected,  so  the  assumptions  made  above  are  well  justified.  It  should  be 
noted  that  in  the  derivation  of  Equation  (5.1c),  the  term  g Qg  was  added  and  subtracted  to 
obtain  the  given  form  of  the  source  term.  As  a result  of  this,  the  pressure  shown  in  the 
momentum  equations  explicitly  contains  the  hydrostatic  component,  thus, 

P = P*  + Qogy  (5.4) 

where  p*  represents  the  pressure  in  the  original  equations  of  Chapter  2,  which  is  due  solely 


(5.1b) 
(5.1c) 
(5. Id) 
(5.1e) 


to  the  fluid  motion. 
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Discretization  of  the  governing  equations  follows  along  the  same  lines  as  that 
detailed  in  Chapter  2,  the  only  additions  being  the  incorporation  of  the  unsteady  terms  in  each 
equation  and  the  buoyancy  term  in  the  v-momentum  equation.  In  this  work,  we  have  adopted 
the  use  of  the  fully  implicit  scheme  (Patankar  1 980)  for  time  marching,  which  is  0(A  t) . For 
the  general  variable  (f> , representing  either  u,  v,  T,  or  S,  the  discrete  form  of  the  governing 
equations  can  be  written  as 

Aj/pp  — A fj>i  + Q (5.5) 

i = E,W,N,S 


where  the  A,-  terms  on  the  right-hand  side  of  the  equation  are  the  same  as  those  shown  in 
Equations  (2.3a-d).  The  term  AP  is  now  given  by 


Ap  = Ap  + Aw  + An  + As  4-  AP  (5.6) 

where  AP°  is  the  contribution  from  the  unsteady  term,  which  takes  the  following  values  for 
the  momentum,  temperature,  and  salinity  equations. 


uorv  momentum : 


o 


A = 

Ap  - 


At 


(5.7a) 


T or  S : 


o AxAy 
p ~ At 


(5.7b) 


For  the  momentum  equations,  the  fluxes  F and  D used  in  defining  the  coefficients  A,  are 
identical  to  those  shown  in  Equations  (2.4)  and  (2.5).  For  the  temperature  and  salinity 
equations,  since  the  density  has  already  been  moved  to  the  right-hand  side  of  the  equation 
to  obtain  the  diffusivity  coefficient,  the  corresponding  forms  of  the  fluxes  are  given  as 

Fe  = ue  Ay  Fw  = uw  Ay  Fn  = vn  Ax  Fs  = vs  Ax  (5.8a) 


a4y  D _ aAy 
{dx)e  w (dx)w 


_ aAx 

(<5y)„ 


, _ aAx 

0y)s 


(5.8b) 


where  a represents  the  appropriate  diffusivity  of  either  temperature  or  salinity.  The 
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generalized  source  term  Q,  shown  in  Equation  (5.5)  takes  the  following  values  for  the 
different  equations 

u:  Q = (Pw  ~ Pe)Ay  + Ap°(pp°  (5.9a) 

v : Q = (ps  ~ Pn)Ax  + Ap°(pp°  + QQglfijiTp  - T0)  + fts(SP  - S0)]AxAy  (5.9b) 

Tor  S:  Q = AP0(pP0  (5.9c) 

In  Equations  (5.9),  the  term  (pP°  represents  the  value  of  the  dependent  variable  from  the 
previous  time  step.  The  terms  TP  and  SP  represent  the  temperature  and  salinity  evaluated 
at  the  center  of  the  v control  volume,  and  are  obtained  by  linear  interpolation  from  the 
staggered  grid  values  above  and  below,  since  T and  S are  stored  in  the  same  physical  location 
as  the  density  and  pressure. 

In  the  following,  we  first  provide  a discussion  of  the  composite  grid  procedure  for 
natural  convection  flows,  where  the  primary  emphasis  is  on  the  development  of  an 
appropriate  interface  treatment  for  the  auxiliary  equation  for  temperature  and  salinity. 
Following  the  discussion  of  the  basic  composite  grid  method,  computations  of  the 
time-dependent,  double-diffusive  flow  in  a solar  pond-type  configuration  are  presented.  The 
main  point  of  interest  to  be  demonstrated  here  is  the  effectiveness  with  which  composite 
grids  can  be  utilized  to  track  localized,  time-evolving  flows.  Finally,  some  conclusions  are 
summarized. 


5.2  Composite  Grid  Method  for  Natural  Convection  Flows 

Since  the  continuity  equation  detailed  in  section  5.1  is  identical  to  that  previously 
used  and  the  momentum  equations  differ  only  by  the  presence  of  additional  source  terms, 
no  changes  are  required  to  the  original  internal  boundary  treatment  discussed  in  Chapter  3. 
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Due  to  the  presence  of  source  terms,  the  original  statement  of  global  conservation,  Equation 
(3.3),  becomes 


where  the  limits  of  i and  j depend  on  the  particular  momentum  equation;  however,  this  does 
not  affect  the  final  condition  for  global  conservation,  Equation  (3.5),  since  the  sources  are 
volumetric  quantities  and  do  not  contribute  to  the  summation  of  the  fluxes  along  the  internal 
boundaries.  With  the  inclusion  of  the  temperature  and  salinity  equations,  we  must  also  now 
ensure  that  the  temperature  and  salinity  fluxes  are  conserved  across  the  internal  boundaries. 
Since  the  temperature  and  salinity  variables  are  stored  in  the  same  physical  locations  as  the 
pressure  on  the  staggered  grid,  the  global  conservation  boundaries  for  both  equations 
correspond  to  the  global  conservation  boundary  for  the  continuity  equation  (which  is  the 
physical  grid  boundary),  and  thus,  a similar  procedure  for  explicitly  conserving  the 
temperature  and  salinity  fluxes  across  the  internal  boundaries  can  be  easily  devised. 
However,  with  the  original  explicit  conservation  procedure,  whereby  the  flux  boundary 
conditions  for  the  temperature  and  salinity  control  volumes  along  the  internal  boundaries  are 
computed  entirely  from  the  neighboring  block,  an  arbitrary  jump  is  allowed  at  the  boundary 
which  cannot  be  detected.  This  non-physical  jump  at  the  interface  between  blocks  can  be 
clearly  seen  from  the  multiblock  solution  of  an  equivalent  one-dimensional  model  problem. 

Consider  the  solution  of  the  convection  diffusion  equation 


(5.10) 


pdT  = d^T 


dx  dx2 


(5.11) 


on  the  two-block  domain  shown  in  Figure  5.1,  where  P is  the  Peclet  number.  The  left  and 
right  boundaries  are  held  at  fixed  temperatures  of  1 and  0 respectively,  while  the  temperature 
at  the  center  of  the  domain  is  taken  to  be  TL  for  the  left  block  (block  1)  and  TR  for  the  right 
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block  (block  2),  indicating  a temperature  jump  of  AT  = TL  — TR  . The  exact  solutions 
within  the  two  blocks  can  be  obtained  as 


Tx{x) 


(5.12a) 


T2(x) 


(5.12b) 


Upon  equalization  of  the  total  temperature  fluxes  at  the  center  of  the  domain, 
computed  independently  from  each  of  the  blocks,  i.e., 


(PT-f), 


L.i 


(5.13) 


JC=1 


we  arrive  at  the  following  relation 


= ,-r  (5.14) 

which  states  that  the  total  fluxes  can  be  balanced  while  still  allowing  a temperature  jump 
AT  = Tl  — Tr,  which  is  a function  of  Peclet  number  and  the  values  of  either  TL  or  T R. 
Thus,  for  an  internal  boundary  treatment  based  on  the  explicit  conservation  of  the  total  flux 
from  the  neighboring  blocks,  a temperature  jump,  the  magnitude  of  which  is  dependent  on 
the  local  flow  conditions,  is  evidently  permitted.  It  should  be  noted  that  for  the  isothermal 
flow  cases  considered  previously,  mass  flux  conservation  uniquely  determines  the  normal 
velocity  component  at  the  interface.  Since  the  normal  and  tangential  velocity  components 
at  the  interface  are  coupled  through  the  continuity  equation,  no  jump  is  permitted  for  the 
tangential  velocity  component  from  the  solution  of  the  tangential  momentum  equation,  even 
though  a locally  conservative  flux  condition  is  employed. 

In  order  to  alleviate  this  temperature  jump  problem,  a new  interface  scheme  based 
upon  a linear  interpolation  procedure  with  a globally  conservative  correction  has  been 
implemented.  By  linearly  interpolating  within  the  neighboring  block  for  the  dependent 
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Figure  5.1.  Two-block  domain  for  one- dimensional  convection-diffusion 
equation,  with  temperature  jump  at  interface. 


variable  (T  or  S)  and  using  this  value  in  conjunction  with  values  in  the  block  of  interest  to 
compute  the  total  flux  boundary  condition,  the  proper  continuity  of  the  temperature  and 
salinity  fields  can  be  maintained.  The  global  correction  is  used  to  adjust  the  fluxes  for  each 
of  the  control  volumes  at  the  boundary  so  that  the  total  temperature  and  salinity  fluxes  along 
the  entire  internal  boundary  separating  the  blocks  are  conserved.  To  illustrate  the 
implementation  procedure  for  the  new  interface  scheme,  consider  a typical  temperature 
control  volume  at  a vertical  internal  boundary,  as  shown  in  Figure  5.2.  With  the  original 
interface  treatment,  the  temperature  flux  through  the  right  control  volume  face  is  obtained 
using  information  exclusively  from  the  neighboring  block  to  the  right  via  an  explicit 
conservation  procedure  as  detailed  for  the  conservation  of  mass  in  Section  3.4.  Using  this 
procedure,  an  arbitrary  temperature  jump  cannot  be  detected  since  the  gradient  term  dT/dx 
is  unaffected  by  a constant  shift  in  the  temperature  field  of  the  neighboring  block.  In  order 
for  the  solution  procedure  to  be  able  to  detect  the  jump  and  adjust  accordingly,  the  gradient 
term  must  be  computed  using  information  from  the  blocks  on  both  sides  of  the  boundary. 
Since  the  temperature  at  the  control  point  on  the  left  side  of  the  interface,  TL,  is  already 
known,  we  must  only  obtain  the  temperature  value  TL  via  linear  interpolation  from  the 
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neighboring  block.  With  the  specified  boundary  velocity  value,  uB  (already  obtained  from 
conservation  of  mass),  the  total  temperature  flux  through  the  right  control  volume  face  can 
be  computed  as 


Again,  we  have  chosen  to  locate  the  imaginary  point  TR  a distance  Ax/2  from  the  boundary, 
in  order  that  the  current  procedure  will  degenerate  appropriately  for  a case  with  a perfect 
overlap  and  continuous  grid  lines  across  the  interface.  Since  this  new  procedure  is  no  longer 
locally  conservative,  we  can  no  longer  expect  global  conservation  to  be  automatically 
enforced,  thus  the  need  for  a global  correction.  This  is  implemented  by  computing  the  total 
temperature  flux  through  line  A independently  from  both  blocks,  using  Equation  (5.15)  for 
the  fluxes  out  of  the  right  block  and  then  adjusting  the  temperature  fluxes  for  each  of  the 
control  volumes  along  the  boundary  by  the  required  constant  to  obtain  global  conservation. 

In  the  above  procedure  for  computing  the  fluxes  for  each  temperature  control  volume 
along  the  interface,  the  use  of  linear  interpolation  is  not  essential,  as  any  order  of 


(5.15) 


A 


Ax  | Ax 

2 I 2 


Figure  5.2.  Temperature  control  volume  located  at  a vertical  internal  boundary. 
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interpolation  may  be  used  for  interpolating  within  the  neighboring  block  (even  piecewise 
constant).  More  importantly,  whatever  form  of  interpolation  is  used,  the  total  flux  (strictly 
only  the  diffusion  term)  must  be  computed  using  information  from  the  blocks  on  both  sides 
of  the  internal  boundary  in  order  to  prevent  temperature  or  salinity  jumps  from  occurring. 
In  this  regard,  no  locally  conservative  boundary  treatment  based  on  piecewise  constant 
interpolation  (as  done  in  the  previous  chapters  for  mass  and  tangential  momentum)  can  be 
successfully  applied  (for  the  T and  S equations)  for  multiblock  problems  requiring  the 
solution  of  the  temperature  and/or  salinity  fields,  since  a locally  conservative  treatment  by 
definition  must  compute  the  fluxes  entirely  from  neighboring  block  information. 

To  test  the  new  interface  treatment,  the  natural  convection  flow  in  the  square  cavity 
shown  in  Figure  5.3  is  computed.  The  boundary  conditions  for  the  flow  and  the  grid  block 
arrangement  are  shown  in  Figure  5.3(a).  The  grid  spacing  in  the  horizontal  direction  is  the 
same  for  both  blocks,  while  in  the  vertical  direction  the  grid  spacing  in  block  one  is  half  that 
of  block  two.  For  this  case,  the  overlap  region  consists  of  a single  column  of  cells.  The  flow 
has  been  computed  for  a Rayleigh  number  of  105  based  on  the  height  of  the  cavity.  Figures 
5.3(b)  and  5.3(c)  show  the  vertical  temperature  profiles  in  the  overlap  region  near  the  center 
of  the  cavity  for  the  solutions  based  on  the  original  conservation  treatment  and  the  new 
treatment.  From  these  plots,  it  is  clear  that  the  original  treatment  of  the  total  temperature  flux 
is  unsatisfactory,  allowing  a large  temperature  jump  in  certain  regions,  while  the  new 
treatment  maintains  the  proper  continuity  along  the  entire  interface. 

5.3  Application  to  Double-Diffusive  Flow  in  a Solar  Pond 
5.3.1  Background 

Double-diffusive  convection  encompasses  the  spectrum  of  buoyancy-induced  fluid 
motions  which  occur  when  two  or  more  components  having  different  molecular  diffusivities 


98 


f = 0 

dy 
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Temperature 
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Figure  5.3.  Two-block  simulation  of  natural  convection  in  a square  cavity,  (a) 
Physical  domain  and  boundary  conditions;  (b)  T profiles  along 
vertical  centerline  for  both  blocks  with  local  flux  conservation 
interface  treatment;  (c)  T profiles  with  interface  treatment  based  on 
linear  interpolation  with  global  flux  correction. 


and  making  opposite  contributions  to  the  fluid  density  gradients  exist  simultaneously.  Since 
the  molecular  diffusivities  of  the  components  can  often  differ  by  an  order  of  magnitude  or 
more  (such  as  when  temperature  and  a salt  are  the  two  components  of  interest),  the  motions 
encountered  in  double-diffusive  systems  can  be  quite  varied  and  complex,  even  for  simple 
physical  systems,  as  summarized  by  Shyy  (1994).  No  attempt  will  be  made  here  to  provide 
a complete  background  of  the  nature  of  double-diffusive  flows,  as  many  reviews  can  be 
found  in  the  literature,  including  Turner  (1974,  1979, 1985)  and  Huppert  & Turner  (1981). 
Instead,  we  concentrate  here  on  one  particular  class  of  double-diffusive  flow  in  which  an 
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initially  quiescent  fluid  with  linear  solute  and/or  temperature  stratifications  is  heated  from 
a sidewall.  Interest  in  this  particular  flow  has  stemmed  primarily  from  observations  of  stably 
stratified  convective  layers  in  natural  reservoirs,  such  as  lakes  and  oceans  (Hoare  1966, 
Huppert  & Turner  1980),  as  well  as  man-made  systems,  including  solar  ponds  (Akbarzadeh 
& Manins  1988,  Sherman  & Imberger  1991).  With  regard  to  solar  ponds,  the  appearance  of 
convective  layers  has  been  found  to  substantially  reduce  the  temperature  difference  which 
can  be  maintained  across  the  pond  depth,  and  thus,  the  efficiency  as  a solar  power  collector, 
although  the  influence  of  sidewall  heating  in  the  generation  of  these  layers  is  a subject  of 
debate  (Sherman  & Imberger  1991). 

The  first  experimental  investigations  of  the  lateral  heating  of  a stably-stratified  fluid 
were  performed  by  Thorpe  et  al.  (1969)  and  Chen  et  al.  (1971).  Both  studies  considered  a 
fluid  of  constant  temperature  with  a vertical  salinity  gradient  subjected  to  heating  from  a 
sidewall  held  at  a fixed  temperature  above  that  of  the  initial  temperature  of  the  interior  fluid. 
These  investigations  clarified  the  nature  of  the  flow  at  the  sidewall  as  a stability  phenomenon 
by  showing  that  under  certain  supercritical  heating  conditions,  convection  cells  appear 
simultaneously  along  the  entire  heated  portion  of  the  wall.  In  addition,  the  nature  of  the 
formation  and  propagation  of  the  intrusions  was  observed,  where  it  was  noted  that  the  initial 
convection  cells  merged  to  form  larger  intrusions  as  they  propagated  laterally  into  the  bulk 
of  the  domain.  Chen  et  al.  (197 1)  proposed  a vertical  length  scale  for  the  intrusions,  h,  given 
by 


_ fiT 

Ps  (dS/dy)i 


(5.16) 


which  characterizes  the  height  to  which  a heated  fluid  element  at  the  wall  would  rise  in  the 
initial  density  gradient.  In  the  above  expression  /3T  represents  the  coefficient  of  thermal 
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expansion,  fis  represents  the  coefficient  of  expansion  due  to  salinity,  A T is  the  difference 
between  the  heated  sidewall  temperature  and  the  initial  temperature  of  the  interior  fluid,  and 
(dS/dy)i  is  the  initial  vertical  solute  gradient,  as  designated  by  the  subscript  i.  Based  on  this 
length  scale,  a critical  Rayleigh  number,  defined  by 


Rac 


gPT  ATh3 
vaT 


(5.17) 


of  15000  ±2500  was  established,  below  which  lateral  intrusions  would  not  form.  In 
Equation  (5.17),  v and  aT  are  respectively,  the  kinematic  viscosity  and  thermal  diffusivity 
of  the  working  fluid. 

Early  investigations  involving  sidewall  heating  focused  on  the  formation  of 
intrusions  due  to  heating  from  a constant  temperature  wall;  however,  in  many  circumstances, 
a constant  lateral  heat  flux  at  the  wall  is  the  appropriate  boundary  condition.  By  performing 
a nondimensional  analysis  of  the  governing  equations  and  assuming  that  the  length  scales 
in  both  the  horizontal  and  vertical  directions  are  identical  at  initiation,  Narusawa  and 
Suzukawa  (1981)  showed  that  in  addition  to  the  Prandtl  number  (77j  = v/aT)  and 
diffusivity  ratio  (Lewis  number  = IJ2  = aT/as  ),  the  third  nondimensional  parameter  is 
given  by  the  following: 


, _ ~ Pjjq/k) 

3 Ps{dS/dy)i 


(5.18) 


where  q represents  the  constant  applied  lateral  heat  flux  at  the  wall,  k is  the  thermal 
conductivity  of  the  fluid  and  the  other  quantities  are  as  previously  defined  in  Equation  (5.16). 
Since  they  considered  a constant  temperature  fluid  with  only  an  initial  linear  solute 
stratification,  this  parameter  can  be  physically  interpreted  as  the  ratio  of  the  horizontal 
density  gradient  at  the  vertical  wall  to  the  initial  vertical  density  gradient.  For  a constant 
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sidewall  heat  flux,  this  third  parameter  was  clearly  identified  as  the  appropriate  stability 
parameter;  however,  the  critical  value  at  which  intrusions  would  develop  was  found  to  be 
dependent  on  the  diffusivity  ratio.  Three  different  solutes  were  tested,  namely,  CUSO4 , NaCl 
and  HC1,  for  which  the  critical  values  of  the  stability  parameter  were  found  to  be  0. 13, 0.28, 
and  0.76,  respectively,  indicating  that  the  critical  value  varies  approximately  linearly  with 
respect  to  the  diffusivity  ratio.  Unlike  the  case  of  a constant  temperature  sidewall  heating, 
no  natural  length  scale  for  the  height  of  the  initial  intrusions  presents  itself  for  a constant  heat 
flux  boundary  condition.  The  length  scale  used  by  Narusawa  and  Suzukawa  (198 1)  in  their 
nondimensional  analysis  was  found  by  them  to  be  an  inappropriate  measure  of  the  height  of 
the  initial  intrusions.  For  values  of  773  between  0.13  and  30.0  the  value  of  the  height  of  the 
initial  intrusions  nondimensionalized  by  this  length  scale  was  found  to  steadily  increase  from 
approximately  7 to  350. 

An  experimental  study  of  double-diffusive  systems  containing  initial  linear  vertical 
stratifications  of  both  temperature  and  salinity  was  recently  carried  out  by  Schladow  et  al. 
(1992).  In  their  experiments,  they  used  two  nondimensional  quantities  to  characterize  the 
flows,  one  the  Rayleigh  (or  buoyancy)  ratio,  given  by 


, = P&S/dy)t 
e /3j(dT/dy)i 


(5.19) 


which  is  formed  from  the  vertical  saline  and  thermal  Rayleigh  numbers  defined  by  the  initial 
salinity  and  temperature  gradients  (again  denoted  by  the  subscript  i)  which  are  both  negative, 
and  the  other  a lateral  ratio,  defined  by 


_ ~ fiiiq/k) 

1 (-PjdT/dy  +/3sdS/dy)i 


(5.20) 


The  Rayleigh  ratio,  Re , provides  a measure  of  the  gravitational  stability  of  the  system,  larger 
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values  indicating  a higher  gravitational  stability.  The  lateral  ratio,  R , , was  proposed  as  a 
modification  to  the  nondimensional  parameter  IJ3  used  by  Narusawa  and  Suzukawa  (1981) 
which  explicitly  takes  into  account  the  two-component  vertical  stratification.  The  quantity 
R j provides  a measure  of  the  forcing  at  the  heated  sidewall  and  can  also  be  interpreted  as 
a ratio  of  the  lateral  density  gradient  at  the  sidewall  to  the  initial  vertical  density  gradient. 
Based  on  these  two  parameters,  a series  of  experiments  was  conducted  using  a NaCl-based 
salt  solution,  aimed  at  delineating  the  physical  characteristics  of  the  intrusions  which 
developed  under  various  states  of  gravitational  and  lateral  stability. 

Three  classes  of  intrusions,  denoted  as  I,  II  and  III,  were  clearly  identified.  Class  I 
intrusions,  which  developed  under  conditions  of  relatively  high  gravitational  and  lateral 
stability  (7?^,  >5  , Rx  ~ 1)  were  characterized  by  relatively  quiescent,  well-defined  cells 
less  than  10mm  thick.  The  internal  structure  of  these  intrusions  consisted  of  a highly  stable 
vertical  temperature  stratification  and  a well-mixed  salinity  field.  The  intrusions  propagated 
at  speeds  less  than  10  cm/h  during  the  period  in  which  the  heat  flux  at  the  wall  was  applied 
and  following  removal  of  the  heat  flux,  they  were  observed  to  quickly  diffuse  away.  Similar 
characteristics  were  observed  by  Narusawa  and  Suzukawa  (198 1)  for  cases  where  the  lateral 
stability  parameter  773  < 1 . In  their  case  RQ  = °°  since  they  considered  a constant 
temperature  fluid.  Intrusions  categorized  in  class  II  developed  under  conditions  of  relatively 
low  gravitational  and  lateral  stability  (2.3  < < 5 , R\  > 1).  These  intrusions  were 

observed  to  be  larger  ( ~ 10-40  mm),  more  dynamic,  and  less-defined  than  class  I 
intrusions.  In  contrast  to  class  I intrusions,  the  vertical  salinity  stratification  was  found  to  be 
slightly  unstable;  however,  the  temperature  field  was  still  highly  stable.  In  addition, 
intrusions  categorized  in  class  II  propagated  at  much  higher  speeds  ( ~ 10  - 30  cm/h)  than 
class  I intrusions;  however,  like  class  I intrusions,  they  also  diffused  away  immediately 
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following  removal  of  the  sidewall  heating.  Class  III  intrusions  developed  under  the 
conditions  of  lowest  gravitational  and  lateral  stability  (Re  ~ 2 , > 1)  and  were  found 

to  be  the  most  dynamic  of  the  three  classes,  propagating  at  speeds  greater  than  25  cm/h.  Due 
to  stronger  convective  motions  within  these  intrusions,  both  the  temperature  and  salinity 
fields  were  found  to  be  well-mixed. 

While  the  physical  appearances  of  class  II  and  class  III  intrusions  were  observed  to 
be  very  similar,  the  defining  characteristic  of  class  III  intrusions  was  their  ability  to  continue 
propagating  long  after  the  removal  of  the  sidewall  heating.  In  this  situation,  sidewall  heating 
serves  as  a triggering  mechanism  for  the  instability,  but  after  removal  of  the  sidewall  heat 
flux,  the  intrusions  are  somehow  capable  of  tapping  off  of  the  potential  energy  stored  in  the 
destabilizing  component  (temperature  in  this  case).  This  observation  of  self-propagation 
could  have  some  implications  regarding  the  development  and  persistence  of  convective 
intrusions  in  poorly-maintained  solar  ponds;  however,  in  most  solar  ponds  the  gravitational 
stability  parameter  Rg  is  held  at  such  high  values  (>  5)  that  class  III  intrusions  are  unable 
to  form. 

Regarding  the  merging  process,  Schladow  et  al.  (1992)  noted  that  the  physical 
mechanisms  involved  in  the  merging  process  for  class  I and  class  II  intrusions  appear  to  be 
fundamentally  different.  Following  the  formation  of  the  initial  intrusions  for  a class  I flow, 
merging  seems  to  take  place  via  a process  whereby  weaker  intrusions  are  forced  back  to  the 
heated  wall  by  a blockage  effect  from  the  intrusions  above  and  below,  with  no  initial 
exchange  of  fluid  between  neighboring  intrusions.  Conversely,  the  merging  process  for  class 
II  and  class  III  flows  occurs  via  specific  events,  in  which  jets  of  fluid  near  the  heated  wall, 
originating  from  some  intrusions,  penetrate  into  neighboring  intrusions  from  below,  causing 
a breakdown  of  the  interface  separating  the  intrusions  and  leading  to  merger.  Similar 
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observations  of  the  merging  processes  discussed  above  have  also  been  made  by  Tanny  and 
Tsinober  (1988). 

Detailed  numerical  simulations  of  the  development  and  propagation  of  intrusions 
due  to  sidewall  heating  have  not  appeared  in  the  literature  to  date,  most  likely  due  to  the 
excessive  computational  resources  required  to  capture  these  flows.  The  physical  time  scales 
of  the  fluid  motion  are  usually  0(0.  Is),  while  the  time  required  to  achieve  full  formation  of 
the  intrusions  and  to  observe  the  merging  processes  is  0(10  min)  or  more  for  slowly 
propagating  intrusions.  Combined  with  the  high  grid  resolutions  required  to  accurately 
capture  the  sharp  velocity,  temperature,  and  salinity  gradients  located  in  the  interfaces 
between  the  intrusions,  the  problem  is  a formidable  one. 

In  addition  to  their  experimental  investigation,  Schladow  et  al  (1992)  conducted  a 
single  numerical  simulation,  solving  the  unsteady  mass,  momentum,  energy,  and  species 
conservation  equations;  however,  they  were  primarily  concerned  with  the  initiation  process 
and  thus  the  computations  were  not  taken  to  the  point  where  full  development  of  the 
intrusions  could  be  observed.  Their  physical  domain  consisted  of  and  8 cm  long  by  2.91  cm 
high  portion  near  the  heated  wall,  which  was  spanned  by  a 98x98  computational  mesh,  with 
free-slip  horizontal  boundaries  and  no-slip  vertical  boundaries.  The  initial  vertical 
temperature  and  salinity  profiles  were  taken  to  be  linear  and  the  horizontal  boundaries  were 
held  at  the  constant  flux  values  required  to  sustain  the  initial  gradients.  With  the  constant  heat 
flux  boundary  condition  imposed  at  the  heated  wall,  their  numerical  simulation  corresponds 
to  a situation  where  Re  = 5.0  and  Rx  = 0.5  . Unfortunately,  this  numerical  simulation  did 
not  correspond  to  any  of  their  experimental  cases,  and  thus,  direct  comparison  between  the 
two  could  not  be  made. 

Here,  direct  simulations  of  two  of  the  experimental  cases  run  by  Schladow  et  al. 
(1992)  are  performed  using  the  composite  grid  method  which  has  been  developed.  At  the 
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current  time,  it  is  still  not  feasible  to  simulate  the  flow  in  the  entire  experimental  apparatus; 
however,  this  is  not  required.  By  direct  simulation,  it  is  implied  here  that  systems  with 
double-diffusive  fronts  composed  of  multiple  convection  cells  are  considered.  Although  the 
overall  physical  sizes  of  the  configurations  computed  here  are  somewhat  smaller  (by  a factor 
of  two  in  the  vertical  direction,  and  a factor  of  twenty  in  the  horizontal  direction),  the 
configurations  are  large  enough  to  simulate  flows  with  intrusion  fronts  composed  of  a 
sufficient  number  of  convection  cells  in  isolation  from  the  influence  of  the  physical 
boundaries,  that  direct  comparisons  with  the  experimental  results  can  be  made  regarding 
keys  points  of  interest  including  the  intrusion  front  propagation  speed,  the  internal 
characteristics  of  the  intrusions  (primarily  temperature  and  salinity  profiles),  and  details  of 
the  merging  process. 

5.3.2  Grid  Refinement  Strategy 

Due  to  the  localized  nature  of  the  growth  and  development  of  the  intrusions  at  the 
heated  sidewall,  use  of  the  composite  grid  method  with  a mesh  refinement  strategy  based  on 
the  local  addition  and  deletion  of  grid  lines  provides  an  efficient  and  accurate  solution 
procedure  (for  other  conditions  fixed,  such  as  the  convection  scheme,  etc.).  Figure  5.4  shows 
a typical  time  sequence  for  the  development  of  class  I intrusions  over  a portion  of  the  heated 
wall.  At  initiation,  convection  cells  appear  nearly  simultaneously  at  the  wall,  and  quickly 
merge  to  form  the  initial  intrusions.  The  intrusions  are  characterized  by  a relatively  high 
speed  clockwise  rotation,  in  which  heated  fluid  flows  away  from  the  wall  along  the  top  of 
the  intrusions,  cooling  as  it  moves  along,  and  then  returns  along  the  bottom  of  the  intrusions. 
This  process  is  evidenced  in  the  characteristic  downward  sloping  appearance  of  the 
intrusions.  Typical  fluid  speeds  are  0(1  mm/s)  while  the  propagation  speed  of  the  front  is 
only  0(10  cm/h).  Since  the  intrusions  are  abutting,  sharp  velocity  gradients  exist  along  the 
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interfaces  between  the  intrusions.  In  addition,  the  intrusions  are  usually  well-mixed  in  terms 
of  temperature  and  salinity  and  are  separated  at  the  interfaces  by  sharp  temperature  and 
salinity  gradients.  In  regions  away  from  the  intrusions,  flow  velocities  are  much  smaller  and 
the  temperature  and  salinity  profiles  remain  very  close  to  the  initial  stratifications. 

From  the  physical  description  given  above,  it  is  clear  that  a very  fine  grid  resolution 
is  required  in  the  region  near  the  intrusions,  while  a much  coarser  grid  is  sufficient  away  from 
the  intrusions.  A single  grid,  clustered  in  the  region  near  the  intrusions  might  be  sufficient; 
however,  due  to  the  preferred  directionality  of  the  flow,  in  which  a very  fine  vertical  spacing 
must  be  used  to  capture  the  nearly  horizontal  interfaces  between  the  intrusions,  this  choice 
is  not  optimal,  since  the  fine  vertical  resolution  in  the  near  wall  region  will  also  be  imposed 
on  the  far-field  flow.  An  adaptive  grid  strategy  using  a single  grid  within  a curvilinear 
coordinate  framework  could  be  used  to  optimize  the  vertical  resolution  in  the  near-wall 
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region;  however,  many  grid  points  would  still  be  wasted  in  the  bulk  of  the  domain,  since  the 
vertical  resolution  requirements  in  these  two  regions  are  truly  disparate. 

In  our  simulations,  a two-block  composite  grid  is  used  to  track  the  evolution  of  the 
sidewall  intrusions.  One  block,  with  a very  fine  horizontal  and  vertical  resolution  is  placed 
in  the  near-wall  region,  extending  from  the  wall  to  a location  just  ahead  of  the  moving 
intrusion  front.  A second  block  with  a relatively  coarse  horizontal  and  vertical  resolution  is 
used  in  the  bulk  of  the  domain.  These  two  blocks  share  a vertical  overlap  region  one  coarse 
cell  thick  (which  corresponds  to  several  fine-grid  cells).  As  the  intrusion  front  moves  into 
the  bulk  fluid,  grid  lines  are  added  to  the  fine  grid  and  removed  from  the  coarse  grid,  so  that 
the  find-grid/coarse-grid  interface  remains  ahead  of  the  intrusion  front.  A visual 
representation  of  this  process  for  the  intrusion  development  sequence  shown  in  Figure  5.4 
is  displayed  in  Figure  5.5.  By  keeping  the  fine-grid/coarse-grid  interface  sufficiently  far 
ahead  of  the  intrusion  front,  we  can  ensure  that  the  inter-grid  transfer  of  information  required 
for  providing  starting  values  for  the  newly  introduced  fine  grid  nodes  can  be  done  in  an 
accurate  manner.  Since  coarse-grid  lines  are  being  deleted,  no  new  information  is  required 
for  the  coarse  grid  upon  re-gridding  near  the  front.  In  the  following  section,  we  give  a 
detailed  presentation  of  the  simulations  which  have  been  performed. 

5.3.3  Numerical  Simulation  Results 

5.3.3. 1 Basics.  Two  numerical  simulations  were  performed,  one  corresponding  to 
the  development  of  class  I intrusions  and  the  other  to  class  II  intrusions.  For  both  simulations, 
a 20  cm  high  by  20  cm  wide  domain  was  used.  In  the  experiments  performed  by  Schladow 
et  al.  (1992),  the  tank  used  was  50  cm  high  and  400  cm  wide.  This  height  is  sufficient  to 
accommodate  approximately  30  class  I intrusions  and  10  class  II  intrusions  along  the  heated 
left  wall,  based  on  the  final  intrusion  thickness  after  completion  of  the  merging  process. 
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Similarly,  a domain  of  20  cm  in  height  should  allow  approximately  12  class  I intrusions  and 
5 class  II  intrusions  to  develop  along  the  heated  wall,  which  should  be  sufficient  to  allow 
comparison  with  the  experimental  results  for  the  key  points  mentioned  previously. 
Concerning  the  apparent  disparity  in  the  tank  widths,  it  is  noted  that  for  the  duration  of  the 
experiments  (completely  through  the  final  merger  process),  the  intrusion  fronts  penetrated 
no  more  than  about  20  cm  into  the  tank,  and  thus,  a large  portion  of  the  tank  remained 
essentially  quiescent,  except  for  the  small  disturbances  created  by  the  presence  of  the 
intrusions.  Since  the  numerical  simulations  performed  here  will  be  terminated  when  the 
intrusion  fronts  are  in  the  approximate  vicinity  of  the  vertical  centerline  of  the  domain,  the 
right  solid  boundary  should  not  significantly  influence  the  motion  of  the  intrusions 
developing  from  the  left  heated  wall. 


Figure  5.5.  Composite  grid  tracking  of  sidewall  intrusions. 
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The  initial  stratifications  of  temperature  and  salinity  were  prescribed  as  follows.  At 
the  top  of  the  domain,  a 3 cm  unstratified  zone  was  set,  followed  by  a linear  stratification 
region  of  14  cm,  and  then  another  3 cm  unstratified  zone  at  the  bottom.  In  the  experimental 
simulations,  unstratified  layers  were  also  employed  at  the  top  and  bottom  of  the  tank  to 
provide  a nearly  constant  zero  flux  condition  for  temperature  and  salinity  for  the  duration 
of  the  experiment.  Since  these  flux  conditions  can  be  exactly  enforced  in  the  numerical 
simulations,  the  unstratified  zones  used  here  serve  only  as  a buffer  region,  effectively 
isolating  the  intrusion  front  from  the  effects  of  the  upper  and  lower  walls.  A schematic  of 
the  domain  used  for  the  numerical  simulations,  showing  the  boundary  conditions  and  initial 
temperature  and  salinity  stratifications  is  shown  in  Figure  5.6.  The  prescribed  sidewall 
heating  fluxes,  specified  initial  temperature  and  salinity  stratifications  and  the 
corresponding  Rayleigh  ratio,  lateral  stability  ratio  and  intrusion  classification  for  the 
numerical  simulations  are  summarized  in  Table  5.1.  These  simulations  correspond  exactly 
(in  terms  of  the  values  of  Rg  and  RJ  to  the  two  experimental  cases  reported  in  detail  in 
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Figure  5.6.  Schematic  of  physical  domain,  boundary  conditions,  and  initial 
T and  S profiles  for  numerical  simulations. 
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Schladow  et  al  (1992).  As  the  mean  stratification  temperature  and  salinity  at  the  start  of  the 
experimental  simulations  were  unknown,  the  numerical  simulations  were  assumed  to  occur 
at  a starting  mean  stratification  temperature  of  20°C  and  a mean  salinity  stratification  of  30 
kgm-3,  for  which  the  physical  parameters  of  water,  such  as  the  kinematic  viscosity,  thermal 
diffusivity  and  mass  diffusivity,  all  assumed  constant,  were  taken  from  Akbarzadeh  and 
Manins  (1988).  Thus,  small  differences  in  both  the  thermal  and  solutal  Rayleigh  numbers 
may  exist  between  the  experimental  and  numerical  simulations  for  each  case,  although  the 
Rayleigh  ratios,  Re,  are  the  same. 


Table  5.1  Parameter  Specifications  for  Numerical  Simulations 


Intrusion 

Side  Wall 

Initial 

Initial 

Rayleigh 

Lateral 

Class 

Heat  Flux 

dT/dy 

dS/dy 

Ratio 

Ratio 

(W/m2) 

(K/m) 

(kg/m4) 

Rq  (Eq.  5.18) 

Rl  (Eq.  5.19) 

I 

91.67 

-25.37 

-129.58 

8.3 

0.9 

II 

91.61 

-7.01 

-13.80 

3.2 

10.8 

The  grid  used  for  the  simulations,  as  mentioned  previously,  is  composed  of  two 
blocks,  a fine  block  in  the  vicinity  of  the  heated  wall  and  a coarse  block  away  from  the  wall. 
In  the  initial  grid,  the  fine  block  consists  of  131x501  nodes,  while  the  coarse  block  has 
162x251  nodes.  In  the  fine  grid  block,  the  grid  lines  have  been  clustered  toward  the  heated 
wall,  resulting  in  a grid  spacing  ofzlx=0.012cm  (0.06  % of  overall  domain  width)  at  the  wall 
(x=0.0)  and  increasing  linearly  to  Z1jc=0.05  cm  at  x=4.0  cm.  The  vertical  grid  spacing  in  the 
fine  block  is  uniform,  with  Zly=0.04  cm.  In  the  coarse  block,  the  grid  lines  are  uniformly 
spaced  in  both  directions,  resulting  in  grid  spacings  ofzdx=0.1  cm  andzly=0.08  cm.  Since 
the  last  two  vertical  columns  of  the  fine  grid  have  a grid  spacing  ofzlx=0.05  cm,  the  overlap 
region  of  the  fine  and  coarse  grids  consists  of  exactly  two  fine  grid  cells  and  one  coarse  grid 


cell. 


Ill 


As  the  solution  advances  in  time,  eventually  the  intrusion  front  will  approach  the 
fine/coarse  grid  interface  of  the  initial  grid  and  re-meshing  will  be  required.  Re-meshing  will 
be  carried  out  here  when  the  fastest  moving  intrusion  approaches  a distance  of  about  1 cm 
from  the  interface.  Since  the  intrusion  fronts  for  class  I and  class  II  flows  propagate  at 
different  speeds,  the  times  at  which  re-meshing  occurs  will  be  different;  however  when 
re-meshing  is  required,  the  following  procedure  is  used.  Upon  the  need  for  re-meshing,  a 4.0 
cm  extension  is  added  to  the  fine  grid,  while  4.0  cm  of  the  coarse  grid  is  removed,  producing 
the  new  mesh.  The  grid  spacings  in  the  newly  added  portion  of  the  fine  mesh  are  taken  to 
be  uniform,  with  values  ofzljc=0.05  cm  andzly=0.04  cm.  This  process  results  in  a new  mesh 
composed  of  a fine  block  with  2 1 1 x50 1 nodes  and  a coarse  block  with  122x251  nodes.  Values 
for  the  solution  at  the  newly  introduced  fine  grid  nodes  are  taken  from  the  coarse  block  via 
linear  interpolation  before  the  coarse  block  is  restructured.  As  the  solutions  progress,  a 
second  re-meshing  may  also  be  required,  in  which  case  a similar  procedure  as  that  described 
above  for  the  first  re-meshing  is  invoked.  Figure  5.7  shows  the  sequence  of  grids  employed 
in  the  simulations. 

For  both  simulations,  a time  step  At  = 0.2 s was  used.  Since  the  maximum  fluid 
velocities  are  0(  1 mm/s),  this  value  was  heuristically  chosen  so  that  between  successive  time 
steps,  the  fastest  fluid  particles  will  have  traversed  at  most  one  grid  cell.  To  verify  this  choice, 
a single-block  grid  simulation  corresponding  to  the  class  I simulation  described  above  was 
performed  with  time  steps  of  At  = 0.1s  and  At  = 0.2 s.  The  grid  consisted  of  301x301 
uniformly  spaced  nodes.  At  a physical  time  of  2 minutes,  vertical  profiles  of  the  field 
variables  were  taken  at  a location  near  the  heated  wall.  Comparison  of  the  profiles  for  the 
two  time  steps  indicated  that  the  solution  appeared  to  be  time-step  independent  near 
At  = 0.2s . In  this  regard,  it  is  noted  that  the  time  step  employed  by  Schladow  et  al.  (1992) 
in  their  numerical  simulation,  At  = 0.0125s,  seems  excessively  small,  since  the  minimum 
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first  remeshing 


second  remeshing 


Figure  5.7.  Sequence  of  grids  employed  in  numerical  simulations. 


grid  spacing  used  there  was  0.03  cm,  which  is  very  close  to  that  employed  in  the  current 
numerical  simulations.  In  the  following,  we  first  present  results  for  the  class  I simulation, 
followed  by  results  for  the  class  II  simulation. 

5. 3. 3. 2 Results  for  class  I simulation.  A sequence  of  contour  plots  showing  the 
development  of  class  I intrusions  from  the  beginning  of  the  numerical  simulation  to  the 
termination  time  of  36  minutes  is  shown  in  Figure  5.8.  Contours  of  streamfunction, 
temperature  and  salinity  have  been  plotted  at  each  chosen  time.  In  each  of  the  contour  plots, 
only  the  portion  of  the  domain  corresponding  to  the  initial  stratified  region 
(3  cm  < y < 17  cm)  has  been  included.  The  initial  intrusions  form  along  the  heated  wall 
at  a time  of  about  10  minutes,  and  by  about  28  minutes,  have  spanned  the  entire  stratified 
region  along  the  wall.  In  agreement  with  the  experimental  simulation,  the  intrusions  form 
exclusively  from  the  bottom  of  the  domain  to  the  top  (with  the  exception  of  some  small 
circulation  regions  near  the  top).  The  characteristic  downward-sloping  appearance  of  the 
intrusions,  observed  in  all  experiments  for  these  flows,  are  also  observed  here.  From  the 
streamfunction  plots,  it  is  apparent  that  class  I intrusions  are  characterized  by  relatively 
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t = 32  minutes 


t = 36  minutes 


(a) 

Figure  5.8.  Solution  contours  at  various  times  showing  development  of  class  I intrusions, 
(a)  Streamfunction;  (b)  Temperature;  (c)  Salinity. 
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t = 24  minutes 


t = 28  minutes 


(b) 

Figure  5.8  --  continued 
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t = 24  minutes 


t = 28  minutes 


(c) 

Figure  5.8  — continued 


116 


quiescent  motions,  producing  a well-defined  layered  structure  with  distinct  intrusions 
separated  by  thin,  nearly-horizontal  interfaces. 

Vertical  profiles  of  the  temperature  and  salinity  for  the  full  domain  height,  taken  at 
a distance  of  3.0  cm  from  the  heated  wall  (15%  of  width  from  wall),  corresponding  to  the 
times  shown  in  Figure  5.8  are  given  in  Figure  5.9.  From  these  plots,  the  step-like  nature  of 
the  intrusions  is  clearly  evident,  where  it  is  observed  that  that  a highly-stable  temperature 
profile  and  a well-mixed  salinity  profile  have  developed  in  each  of  the  intrusions.  Both  of 
these  characteristics  were  observed  in  the  corresponding  experimental  case.  Results  for  the 
overall  size  of  the  intrusions  and  the  front  propagation  speed  also  show  favorable  agreement 
with  the  values  observed  experimentally.  For  the  numerical  simulation,  the  average  height 
of  the  intrusions  measured  at  the  heated  wall  at  a time  of  36  minutes  is  approximately  8 mm. 
The  corresponding  experimental  values,  taken  from  photographic  images  at  similar  times, 
is  about  10  mm.  Using  the  sequence  of  streamfunction  contour  plots  shown  in  Figure  5.8(a), 
the  intrusion  front  speed  is  estimated  to  be  approximately  7 cm/h.  In  this  computation,  only 
the  ten  intrusions  located  about  the  center  of  the  stratified  region  were  used,  since  the 
development  of  the  intrusions  near  the  upper  and  lower  portions  may  have  been  hindered  by 
the  fairly  strong  recirculation  regions  which  exist  in  the  unstratified  layers  at  the  top  and 
bottom  of  the  domain.  In  this  regard,  the  attempt  to  isolate  the  intrusion  front  from  the  effects 
of  the  solid  walls  by  including  unstratified  layers  above  and  below  may  have  resulted  in  a 
greater  disturbance  to  the  upper  and  lower  intrusions  than  would  have  existed  without  these 
layers.  No  value  for  the  average  speed  of  the  intrusions  for  the  corresponding  experimental 
case  was  given,  but  the  value  computed  above  does  agree  favorably  with  the  speed  of  1 0 cm/h 
or  less  which  was  generally  observed  for  class  I intrusions. 

Regarding  the  merger  process,  some  evidence  of  the  blockage  mechanism  described 
earlier,  and  observed  clearly  in  the  experiments,  can  be  seen  in  the  final  four  frames  of  the 
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Temperature 
t = 1 8 minutes 


Temperature 
t = 28  minutes 


Temperature  Temperature 

t = 32  minutes  t = 36  minutes 

(a) 

Figure  5.9.  Temperature  and  salinity  profiles  taken  at  a distance  of  3.0  cm 

from  the  heated  wall  at  various  times,  (a)  Temperature;  (b)  Salinity. 
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(b) 


Figure  5.9  — continued 
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streamfunction  contour  sequence  shown  in  Figure  5.8(a)  for  the  fourth  intrusion  from  the 
bottom.  From  a time  of  24  minutes  to  36  minutes,  a progressive  weakening  of  this  intrusion 
is  clearly  evident  as  it  becomes  overtaken  by  its  upper  and  lower  neighbors.  Consistent  with 
the  description  of  the  merger  mechanism,  no  exchange  of  fluid  between  neighboring 
intrusions  appears  to  initiate  the  process  of  merger,  since  the  interfaces  between  the  intrusion 
and  its  upper  and  lower  neighbors  are  clearly  seen  to  be  maintained  during  the  entire  process. 

5.3. 3. 3 Results  for  class  II  simulation.  A sequence  of  contour  plots  for  the 
development  of  the  class  II  intrusions  in  increments  of  4 minutes  from  the  beginning  of  the 
simulation  to  the  termination  time  of  20  minutes  is  shown  in  Figure  5.10.  Again,  only  the 
portion  of  the  domain  corresponding  to  the  initial  stratified  region  has  been  included.  The 
initial  intrusions  form  along  the  heated  wall  at  a time  of  about  2 minutes,  eventually  spanning 
the  length  of  the  heated  wall  in  the  stratified  region  by  a time  of  5 minutes.  Unlike  the  class 

I intrusions,  which  formed  exclusively  from  the  bottom  of  the  domain  to  the  top,  the  initial 
formation  of  the  class  II  intrusions  occurs  at  both  the  lower  and  upper  portions  of  the 
stratified  zone.  As  time  elapses,  the  characteristic  structure  of  the  flow  becomes  evident, 
where  we  again  have  fairly  well-mixed  regions  separated  by  distinct  interfaces  in  which 
sharp  gradients  of  temperature  and  salinity  exist.  It  is  apparent  from  the  salinity  contours  that 
class  II  intrusions  are  much  more  dynamic  and  much  less  structured  than  class  I intrusions, 
which  is  in  good  agreement  with  the  experimental  observations. 

Figure  5. 1 1 displays  vertical  profiles  of  the  temperature  and  salinity  (again  taken  3.0 
cm  from  the  heated  wall)  for  the  full  domain  height  corresponding  the  the  times  shown  in 
Figure  5.10.  Again,  the  step-like  nature  of  the  intrusions  is  evident,  where  it  is  observed  that 
a highly  stable  temperature  profile  has  developed  in  each  of  the  intrusions.  In  contrast  to  the 
class  I intrusions,  in  which  the  salinity  was  well-mixed,  the  salinity  distribution  in  the  class 

II  intrusions  is  observed  to  be  slightly  unstable.  These  characteristics  for  the  internal 
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t = 4 minutes 


t = 8 minutes 


t = 20  minutes 

(a) 

Figure  5.10.  Solution  contours  at  various  times  showing  development  of  class  II  intrusions, 
(a)  Streamfunction;  (b)  Temperature;  (c)  Salinity. 
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t = 20  minutes 
(b) 

Figure  5.10.  — continued 
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t = 4 minutes 


t = 8 minutes 


t = 20  minutes 

(c) 


Figure  5.10.  --  continued 
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temperature  and  salinity  profiles  were  also  observed  in  the  corresponding  experimental  case. 
Regarding  the  apparent  “jagged”  behavior  of  the  salinity  profiles  compared  to  the 
temperature  profiles,  it  is  noted  that  since  the  Schmidt  number  is  very  large  (Sc  = 675),  the 
fine  scale  salinity  structure  of  the  intrusions  is  essentially  being  convected  with  diffusion 
having  a role  to  play  only  within  small  length  scales. 

Results  for  the  overall  size  and  front  propagation  speed,  again,  show  favorable 
agreement  with  the  values  observed  for  typical  class  II  intrusions.  Based  on  the  general 
observations  of  the  experiments,  the  final  height  for  typical  class  II  intrusions  was  found  to 
fall  in  the  range  10^10  mm.  From  the  salinity  and  temperature  profile  plots  for  the 
experimental  case  corresponding  to  the  current  numerical  simulation,  taken  3 cm  from  the 
heated  wall  at  a time  of  19  minutes,  the  average  height  of  the  intrusions  is  estimated  to  be 
about  40  mm.  Based  on  a similar  analysis  of  the  temperature  and  salinity  profiles  for  the 
numerical  simulation,  an  average  intrusion  height  of  about  35  mm  is  obtained.  Using  the 
intrusion  salinity  sequence  shown  in  Figure  5. 10(c),  the  intrusion  front  speed  was  computed 
to  be  about  27  cm/h.  In  this  computation,  only  the  two  center  intrusions  were  used.  This  value 
also  agrees  very  well  with  the  range  of  10-30  cm/h  observed  in  the  experiments  for  class  II 
flows.  Again,  no  estimate  of  the  front  propagation  speed  for  the  specific  experimental  case 
corresponding  to  the  numerical  simulation  was  given. 

Finally,  in  terms  of  the  merger  process,  it  is  clear  that,  unlike  class  I flows,  class  II 
flows  are  dominated  by  dynamic  vertical  motions  in  the  near- wall  region.  These  motions  are 
evidenced  in  the  time-history  sequence  by  the  sharp  salinity  gradient  regions  which  are 
pulled  up  along  the  wall  and  are  subsequently  folded  back  into  the  interior  of  the  intrusions 
when  they  reach  their  apex.  Although  details  of  the  initial  merger  process  are  not  captured 
in  the  plots,  a vivid  example  of  the  merger  process  can  be  seen  in  the  three  final  frames  of 
the  time  sequence  for  the  two  top  intrusions.  Here,  vertical  motions  from  the  lower  intrusion 
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Temperature 
t = 8 minutes 


Temperature 
t = 16  minutes 


Temperature 
t = 20  minutes 

(a) 

Figure  5.11.  Temperature  and  salinity  profiles  taken  at  a distance  of  3.0  cm  from  the 
heated  wall  at  various  times,  (a)  Temperature;  (b)  Salinity. 
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Salinity 
t = 20  minutes 


Figure  5.11  — continued 


126 


have  broken  down  the  interface  separating  it  from  the  intrusion  above.  This  is  in  contrast  to 
the  merger  process  observed  in  the  previous  section  for  class  I intrusions,  which  was  initiated 
by  stronger  intrusions  blocking  weaker  ones  from  propagating  and  forcing  them  back  to  the 
heated  wall.  For  both  merger  processes,  fluid  is  eventually  exchanged  during  merger; 
however,  with  class  II  flows  the  kinetic  energy  of  the  fluid  in  the  near- wall  region  serves  as 
the  mechanism  for  merger. 


5.4  Concluding  Remarks 

In  this  chapter  some  issues  related  to  the  extension  of  the  current  method  in  handling 
flows  with  additional  scalar  transport  equations  have  been  addressed.  It  has  been 
demonstrated  that  an  information  transfer  strategy  based  on  local  flux  conservation  for 
transport  equations  such  as  the  energy  equation  and  species  equations  allows  a jump  in  the 
solution  across  block  interfaces  which  is  not  physical.  With  the  same  conservation  strategy, 
such  a jump  is  not  permitted  for  the  tangential  momentum  equation  due  to  the  mass 
continuity  constraint.  A new  interface  scheme  based  on  a linear  interpolation  strategy  with 
a globally  conservative  correction  has  been  employed  to  overcome  this  difficulty.  The 
double-diffusive  flow  in  a rectangular  enclosure  has  been  computed  to  asses  the  capability 
of  using  locally  refined  composite  grids  for  handling  time-dependent  fluid  flows.  Results  for 
two  classes  of  time-dependent  thermohaline  intrusions  have  been  obtained  which  agree  very 
well  with  experimental  predictions  in  terms  of  the  overall  physical  size  and  structure  of  the 
intrusions,  and  the  intrusion  propagation  speeds.  In  addition,  the  physical  characteristics  of 
the  merging  processes  for  the  two  classes  have  also  been  observed  which  agree  qualitatively 
with  the  descriptions  given  by  the  experimental  simulations. 


CHAPTER  6 

SUMMARY  AND  SUGGESTIONS  FOR  FUTURE  WORK 


A pressure-based  composite  grid  methodology  for  the  incompressible  Navier-Stokes 
equations  is  presented.  The  primary  issue  of  importance  is  the  transfer  of  information 
between  the  different  blocks  of  the  composite  grid.  A conservative  interface  strategy  based 
upon  the  local  conservation  of  mass  and  tangential  momentum  flux  has  been  developed 
which  remains  true  to  the  fundamental  nature  of  the  staggered  grid,  in  which  no  artificial 
pressure  boundary  conditions  are  required.  For  simply-connected  domains,  it  has  been 
shown  that  the  explicit  conservation  of  the  normal  momentum  flux  across  internal 
boundaries  provides  the  means  for  eliminating  the  arbitrary  constants  of  the  pressure  field 
within  the  various  blocks  which  arise  due  to  the  use  of  normal  velocity  component  boundary 
conditions  for  the  pressure  correction  equation.  For  multiply-connected  regions,  the  explicit 
balance  of  normal  momentum,  within  the  context  of  a pressure  cycling  strategy,  has  been 
shown  to  play  a more  fundamental  role  in  ensuring  that  unique  solutions  obeying  all  of  the 
conservation  laws  are  obtained. 

In  order  to  provide  a general  capability  for  complex  configurations,  as  well  as  to 
provide  a logical  framework  for  implementing  the  conservative  internal  boundary  scheme, 
an  organizational  procedure  has  been  devised  based  on  the  idea  of  boundary  segments,  which 
is  capable  of  handling  configurations  composed  of  an  arbitrary  number  of  overlapping  grid 
blocks.  Two  demonstrative  configurations  have  been  presented  which  show  that  the  current 
procedure  is  very  well  suited  for  internal  fluid  flows  with  multiple  obstacles,  which  by  their 
nature,  are  best  handled  by  the  composite  grid  method. 
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Several  issues  regarding  the  implementation  of  the  multigrid  technique  within  the 
composite  grid  methodology  have  been  addressed.  Using  a new  grid  coarsening  procedure, 
similar  to  that  used  previously  by  others  for  turbulent  flow  computations,  an  efficient  and 
flexible  composite  grid  method  is  devised.  The  composite  multigrid  method  developed  here 
uses  the  composite  grid  procedure  as  a module  which  operates  within  the  outer  framework 
of  a FMG-FAS  multigrid  algorithm,  whereby  at  each  multigrid  level,  the  appropriate 
composite  grid  problem  is  solved.  Model  problem  testing  for  two  multiply-connected 
configurations  has  demonstrated  that  the  current  composite  multigrid  method  can  maintain 
the  same  level  of  effectiveness  as  its  single-grid  counterpart,  even  for  highly  complex 
configurations. 

Finally,  some  issues  related  to  the  extension  of  the  current  method  in  handling  flows 
with  additional  scalar  transport  equations  have  been  addressed.  It  has  been  demonstrated  that 
an  information  transfer  strategy  based  on  local  flux  conservation  for  transport  equations  such 
as  the  energy  equation  and  species  equations  allows  a jump  in  the  solution  across  block 
interfaces  which  is  not  physical.  With  the  same  conservation  strategy,  such  a jump  is  not 
permitted  for  the  tangential  momentum  equation  due  to  the  mass  continuity  constraint.  A 
new  interface  scheme  based  on  a linear  interpolation  strategy  with  a globally  conservative 
correction  has  been  employed  to  overcome  this  difficulty.  The  double-diffusive  flow  in  a 
rectangular  enclosure  has  been  computed  to  asses  the  capability  of  using  locally  refined 
composite  grids  for  handling  time-dependent  fluid  flows.  Due  to  the  nature  of  the  particular 
flows  chosen,  in  which  isolated  evolving  flow  features  exist  with  sharp  gradients  in  a 
preferred  direction,  adaptive  composite  grids  provide  a distinct  advantage  over  single  grids, 
allowing  the  overall  accuracy  to  be  improved  while  expending  roughly  the  same 
computational  effort.  Results  for  two  classes  of  time-dependent  thermohaline  intrusions 
have  been  obtained  which  agree  very  well  with  experimental  predictions  in  terms  of  the 
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overall  physical  size  and  structure  of  the  intrusions,  as  well  as  the  intrusion  propagation 
speeds.  In  addition,  the  physical  characteristics  of  the  merging  processes  for  the  two  classes 
have  also  been  observed  which  agree  qualitatively  with  the  descriptions  given  by  the 
experimental  simulations. 

In  this  work,  several  fundamental  issues  regarding  the  development  of  a flexible  and 
efficient  composite  grid  method  for  the  pressure-based  SIMPLE  algorithm  have  been 
addressed.  A firm  foundation  has  been  established,  however,  additional  work  is  required 
before  the  current  methodology  can  be  applied  in  a truly  useful  way  to  complex  problems 
of  practical  engineering  interest.  The  current  methodology  has  been  devised  in  the  context 
of  a Cartesian  component  grid  system,  and  is  hence  limiting  in  its  scope  of  application.  In 
order  to  expand  the  capability  to  handle  real  geometries,  a generalized  curvilinear  coordinate 
framework  must  be  adopted.  With  the  introduction  of  curvilinear  coordinates,  the  bulk  of 
the  original  methodology  still  remains  intact,  however,  some  additional  issues  arise  due  to 
the  presence  of  pressure  in  both  the  normal  and  tangential  momentum  fluxes,  as  well  as  the 
appearance  of  cross-derivative  terms  in  the  normal  momentum  flux  term  at  the  internal 
boundaries.  These  issues  are  considered  to  be  minor  ones;  however,  and  do  not  impose  any 
substantial  difficulties  in  extending  the  current  method  to  general  curvilinear  composite  grid 
systems. 

Additional  future  work  should  involve  the  extension  of  the  method  to  compressible 
flows,  as  well  as  the  incorporation  of  some  form  of  turbulence  model  for  computing  high 
Reynolds  number  flows.  For  compressible  flows,  the  only  modification  required  is  to  take 
the  effect  of  pressure  on  density  into  account  in  obtaining  the  correct  form  of  the  pressure 
correction  equation.  As  far  as  the  internal  interface  procedure  is  concerned,  no  modifications 
are  required  to  extend  the  method  to  compressible  flows.  However,  it  will  be  interesting  to 
see  how  the  current  interface  strategy  performs  in  the  presence  of  sharp  gradients  such  as 
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shock  waves  passing  through  internal  boundaries.  In  this  regard,  the  current  methodology 
should  also  be  modified  to  incorporate  high-resolution  such  as  those  detailed  in  Thakur 
(1993).  Regarding  turbulence,  since  the  governing  equations  cannot  be  cast  into  a strong 
conservative  form,  such  as  with  the  two-equation  k-e  model,  and  since  source  terms 
containing  higher-order  mixed  partial  derivatives  are  often  present,  some  investigation  of 
the  method  of  information  transfer  used  at  the  internal  boundaries  for  the  turbulence 
equations  is  also  in  order. 
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